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CHAPTER 1 


INTRODUCTION 


An important class of Shuttle deployed payloads consists of canti- 
levered beam-like structures with massive tip bodies- This report pre- 
sents analytical dynamics models for the planar motion of such Orbiter- 
payload systems (Figure 1.1) . The models are specifically intended for 
use in stability studies of the Orbiter flight control system. Well 
established engineering approximations are invoked m the interests of 
simplicity and tractability . The format is a stepwise progression of 
mechanics problems each providing useful results and insight and forms a 
basis from which to address more complicated situations. 

The payload beam- like structures are taken throughout to be long 
slender uniform beams capable of transverse bending only. A sequence of 
classical beam vibration eigenvalue problems are examined, namely: a 

cantilever with tip mass, a cantilever with tip body and an unconstrained 
beam with rigid bodies at each end. In each case the frequency equa- 
tions, eigenfunctions and orthogonality relations are derived. The anal- 
ytical treatment permits the free vibration characterization in terms of 
a minimum number of dimensionless parameters. As a precursor of the 
ultimate problem, the forced vibration of a cantilevered beam with tip 
body subject to base acceleration is studied. The exact solution to the 
nonhomogeneous partial differential equation with time dependent boundary 
conditions is presented. Natural "modal parameters" are defined and 
important identities in terms of these quantities derived. An approxi- 
mate solution using the assumed modes method proves revealing and serves 
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to check the exact solution. Finally, the planar dynamics of an elastic 
beam with rigid bodies affixed to the ends is addressed. This model of 
the Orbiter-payload system is capable of arbitrary motion in the plane 
accompanied by small elastic deformations. External forces and torques 
acting on the rigid bodies are accommodated. The vehicle motion equa- 
tions are derived for two disparate modal expansions of the beam deflec- 
tions. The equations of motion are cast in a first order form suitable 
for numerical integration. FORTRAN computer programs implementing the 
motion equations are given. 
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CHAPTER 2 


EIGENVALUE PROBLEMS 


2.1 Natural Frequencies and Mode Shapes of a Cantilevered Beam with 
Tip Mass 

The partial differential equation governing the transverse vibra- 
tion of a beam is given by 


-4 .2 

3 U , d u 

^ 4 + ^ 2 

3x at 


0 


This equation assumes a uniform distribution of stiffness El and mass per 
unit length p . The beam is clamped at x = 0 

u(0,t) = (0,t) = 0 for all t > 0 

O X 

at the other end of the beam: x = l, there is a point mass m^ (Fig- 

ure 2.1). 
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Figure 2.1. Clamped beam with tip mass. 
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(Aft) 


0 


a 2 u 


9x 


and 


.3 3 2 

El -2-Si (A,t) - m -2-H (A #t ) 

8x at 

Assume a solution of the form u = <{>{x)e 1Ut . We then obtain the following 
boundary value problem for the mode shapes <j> 

j 4 , 

= 0 «j»(0) = <f>' (0) =0 

dx 

2 m t 

(A = pu /El) <j)"(A) = 0, <j) ' ' ' (A) = - ~ \<p(Z) 

( 2 - 1 ) 

The orthogonality condition m this case can be arrived at by a consider- 
ation of the kinetic energy expression for the system 

o 

1 • 2 i 2 

T = — p / u(x,t) dx + — m u ( A,t) 

0 

Expanding u(x,t) m a series of modes 


u 



(x)q k (t) 


5 



the kinetic energy can be written as a quadratic form m {q } 

K 


co oo / a 

T = j Z Z )p / 

n=l m=l ( 0 


We suspect that the natural modes will diagonalize this form, i.e. 


£ 

P / <p (x) <j> (x) dx + m 4> (£)<£> (£) = 0 (2-2) 

n m t n m 

where <f> (x) and d> (x) are characteristic functions corresponding to dis- 
n m 

tinct characteristic numbers. Indeed this orthogonality condition can be 

verified directly by appealing to the differential equation and boundary 
* 

conditions (2.1) . 

4 

If we set X = a the general solution of the differential equation 


<p (x) = c 1 sm ax + c 2 cos ax + c^ sinh ax + c^ cosh ax 

(X = 0 is not an eigenvalue) 


The boundary conditions at x = 0 require 


d> (x) 4> (x) dx + m. 4> (£) 
n m t n 


(A). 


m 


'q q 

n m 


C 2 + C 4 


0 


and 


Applying the boundary conditions at x = £ we are lead to the frequency 
equation 


* 

This alternate method of arriving at the orthogonality condition was 
first suggested by P.C. Hughes of UTIAS. 
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m t “t 

cos al + cosh al + — alsinh al - sin a l) sinh al - sin al + — a (cosh al - cos al) 

P p 


sin al + sinh al 


cos al + cosh al 


= 0 


Expansion and simplification yields 
m t 

8 (sin 8 cosh 8 - cos 8 sinh 8) = 1 + cos 8 cosh 8 (2-3) 

where we have introduced the dimensionless parameter 8 = al. Note that 
the roots of this equation only depend upon the ratio of the tip mass to 
the mass of the beam. 

The natural frequencies u> k are given by 


, El .2 

“ k JpSL 4 k 


and the corresponding mode shapes ( x) are 

JC 


(2-4) 


+k W 


A k [ (cos 8 k + cosh 8 k )(sin 8 k f - sinh 8 k f) 
- (sin 8 k + sinh 0 k ) ^cos 8 k j - cosh 8 k j)J 


(2-5) 


If we set m = 0, the orthogonality condition (2-2) and the frequency 
equation, Eq. (2-3) , reduce to those for an ordinary clamped-free beam. 

In the limit as -*■ 00 the boundary conditions at x = i approach those of 
a pinned connection: zero displacement and zero moment. And indeed the 

frequency equation, Eq. (2-3) , takes the form appropriate to a canti- 
levered-pmned beam (tan 8 = tanh 8) . 
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Natural Frequencies and Mode Shapes of a Cantilevered Beam with 
Trp Body 

Figure 2.2 depicts the situation m the undeformed position. A 
uniform beam of mass density p, bending stiffness El, and length i coin- 
cides with the x axis. At the tip (P) a rigid body is attached of mass m 
and inertia J (about P) . The distance between P and the rigid body mass 
center is c and this directed line segment makes an angle y with the 
positive x axis. 

v 



Figure 2.2. Clamped beam with tip body - undeformed. 


Figure 2.3 shows the system m a deformed position. Note that the 

tip body is 'rigidly' attached at P. Denote the inertial velocity of P 

by v and the angular velocity of a body frame (fixed in the tip body) by 
P 



Figure 2.3. Clamped beam with tip body - deformed. 
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u(£,t) 3, 


(£,t)k 


v 

P 


0) 

p 


9 2 u 
9 x 3 t 


where u(x,t) is the elastic deflection of the beam and 1, 3, k are the 
unit vectors along axes x, y, z. Denoting the velocity of the mass cen- 
ter of the tip body by v and the vector from p to the mass center by c 

w 

we have 


V + 0 ) x c 
P P 


observing that c = c cos (y + 0 ) 1 + c sin (y + 0^) 3 , where 0^ is the 

angle between the positive x-axis and the beam tip tangent. Noting that 

0 I << 1 so that sin 0 s 0 , cos 0 £ 1 we arrive at 

P 1 P P P 


1 — 


3 2 

t) + C 9x g t - U,t) (cos y - 0 p Sin y) 3 


- c 


a 2 

9x g t (Z,t) (sin y + 0 p cos y)i 


The nonlinear term 0 


2 

d ^ (Z f t) =00 will be dropped 
p 3 x 3 t P P 


- [$?<*• 


t) + c cos y 


3 x 


at 


2 

■ c sin y (Zft)I 


( 2 - 6 ) 


The linear and angular accelerations of the tip body follow directly by 
differentiation of the above expressions. It will be observed that in 
general the mass center of the tip body will have a component of acceler- 
ation along the x axis 
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-c sm y 


U,t) 


This implies a force acting on the beam along the x direction. 

A problem arises at this point if we wish to accommodate the af- 
fect of axial loading on the transverse bending of the beam since non- 
linear equations would result m the context of the present investiga- 
tion. For simplicity we therefore assume that y = 0 throughout the re- 
mainder of this section. The mass center offset of the tip body is 
therefore restricted to be directed axially from p and this results m no 
axial force being applied to the beam. 


= ^-tr(A,t) + c - 3 - — • U,t) ; 

_3t 3x 3t 


b) 

P 


- 3 U 7 U,t)k 

9x at 


(2-7) 


a is the acceleration of the mass center of the tip body and to is the 

w P 

angular acceleration of the tip body. 

In order to write the boundary conditions for u(x,t) at the end- 
point x = £ we consider a free body diagram of the tip body . As indi- 
cated in Figure 2.4 the beam exerts a force S directed along the y axis 
at p and a moment M directed along z. The equation of motion of the tip 
body along the y axis is 
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Figure 2.4. Free body diagram of tip body. 


The shearing force in the beam at x = l is given by S = El 
This gives one of the required boundary conditions 


3 3 u 

3x 3 


a,t) . 


El (A,t) 

3x 


m 


3 2 u 

3t 2 


(£,t) + c 


3 3 u 


(A,t) 


3x 3t 


= 0 


( 2 - 8 ) 


The second boundary condition is obtained by writing rotational motion 
equations for the tip body. Let [I] denote the inertia matrix of the tip 
body about its mass center and h its angular momentum both referred to 
body axes at the mass center. 



[I]oi + U) 
-P “P 


mm 

“P 


Dropping the nonlinear term in it follows that 




, 2 3 3 u 

(J - me ) ; 

3x 3t^ 


U,t) 


where J is the moment of inertia of the tip body about an axis parallel 
to z and passing through P. 

From Figure 2.4, the net moment about the mass center is 

M - c x ? 
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M can be calculated from beam theory in terms of stiffness and curvature. 
Specifically 


3 2 

M = -El -2-5- (A,t)k 
9x 


From the work above we readily obtain 


c x S = cm 


ifu 

9t 2 


(A,t) + c 


a 3 
9 u 

9x 9t z 


(A f t) 


Taking the z component of the rotational motion equation 


_d_ 

dt 


h 


M - c x s 


results m the 2nd boundary condition at x = i 


9 2 9 2 

El — ^ (&,t) + me — 7T (£,t) + J 


9x 


9t 


9 3 u 
9x 9t^ 


(Art) = 


(2-9) 


Since the beam is clamped at x = 0 we have the two geometric boundary 
conditions 


u(0,t) =0 and ^ (0,t) = 0 

dX 


( 2 - 10 ) 


The partial differential equation for free vibration is of course 


a 4 a 2 

9 u 9 u 

EI ~4 + P ~2 = 0 

9x 9t 


( 2 - 11 ) 
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We now proceed to solve the partial differential equation (2-11) subject 


to the geometric boundary conditions (2-10) and the natural boundary 
conditions (2-8) and (2-9). Seeking solutions of the form e 1Wt (f>(x) we 
are led to the eigenvalue problem 


d 4 4> 

dx 4 


- H = o 


( 2 - 12 ) 


>"'(«.) + m — [<j> W + c<f>' (£) ] 
P 


= 0 


(2-13) 


<j>' ' (£) - — [mc<f) (£) + J<j>' (£) ] 


= 0 


(2-14) 


<M0) = 4>'(0) = 0 


(2-15) 


Proceeding as m previous sections, the orthogonality condition can be 
arrived at by considering the kinetic energy, T, of the system. 


1 * 2 
T = — f u(x,t) p dx 

0 


1 2 . 
+ — (J - me ) 




Expanding u(x,t) = £<j> (x)q (t) the orthogonality condition is 

k x X 


z 

p / <f> (x) <j> (x) dx + m<(> (Z) <f> (Z) + J<p ' (Z) p ' (Z) 

0 1 D 13 ID 

+ met^ (Z) <p^ (Z) + <P^(Z)<t>j( Z)] = 0 

(2-16) 
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where <f> (x) , 
values . 


$ (x) 


are eigenfunctions corresponding to distinct eigen- 


Returning to the eigenvalue problem, it can be readily shown that 
A = 0 is not an eigenvalue. The general solution of (2-12) can be writ- 
ten as 


<p (x) = c^ sin ax + c 2 cos ax + c^ sinh ax + c^ cosh ax 

(A = a 4 1 0) 

The boundary conditions (2-15) require c^ = -c and c^ = -c^. Applying 
the boundary conditions (2-13) , (2-14) and eliminating c^ and c^ we ar- 
rive at the following simultaneous equations in c^ and c 2 


(sm 0 - sinh 0) + -^r- c0(cos 0 - cosh 0) — (cos 0 + cosh 0) c 

LP P* p J 1 

+ (cos 0 - cosh 0) — c0(sin 0 + sinh 0) + 4 (sm 0 - sinh 0) c = 0 

|_p pJt p J ^ 

(2-17) 


- — (sin 0 - sinh 0) - —r 0(cos 0 - cosh 0) - — (sm 0 + smh 0) c, 

I P p * 0 2 J 1 

+ I- — (cos 0 - cosh 0) + —77 0(sin 0 + smh 0) - (cos 0 + cosh 0) c„ = 0 

L P p£ e 2 J 2 


(2-18) 


where we have introduced the symbol 0 = al. Setting the determinant of 
the system (2-17) , (2-18) to zero gives the transcendental frequency 
equation 
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m* ( J* 


m^c’ 


4 

) 8(1 - 
2m*c*8 2 


cos 3 cosh 
sin 3 smh 


3) + m*3(cos 3 smh 3 - sin 3 cosh 3) 

3 - J*3 3 (sm 3 cosh 3 + smh 3 cos 3) 

+ 1 + cos 3 cosh 3=0 
(2-19) 


where the dimensionless tip body parameters are defined by 


m* 



c* 



The natural frequencies are given by 



and the corresponding eigenfunctions by 


2 

k 


( 2 - 20 ) 


<0 k (x) = j^m*3 k (cos 3 k - cosh 3 k ) - m*c*3 k (sin B k + smh B k ) 

+ sin B k - smh 8 fc ] • (sm 3 k f - smh 3 k f ) 

+ m*3, (smh 3, - sm 3, ) + m*c*3, 2 (cosh 3, - cos 3, ) 

L k k k k k k 

+ cos B k + cosh 3 k J • |cos 3 k ^ - cosh 3 k 

( 2 - 21 ) 

Note: the eigenfunctions given in (2-21) are not normalized. 

2.3. Natural Frequencies and Mode Shapes of an Unconstrained Beam with 
Two Tip Bodies 

Figure 2.5 depicts the situation m the undeformed state. A uni- 
form beam of mass per unit length p, bending stiffness El, and length £ 
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Figure 2.5. Unconstrained beam with rigid bodies 
attached to each end - undeformed. 


lies on the x axis x <_ Z) . A rigid body of mass m^ and inertia 

about its mass center, is attached at x = 0. The mass center offset is 
c^ directed axially. A second rigid body is attached at x = H with asso- 
ciated parameters m^, and c defined similarly. The motion is planar 
and the beam is capable of bending only. 

The partial differential equation for free vibration is 


El 



+ P 



0 


where u(x,t) is the transverse displacement of the beam from the neutral 
axis. The boundary conditions at x = £ can be obtained directly from the 
analysis in section 2.2. Due to subtle sign changes the kinematics of 
body 1 will be derived here so as to obtain the correct boundary condi- 
tions at x = 0. Denote the inertial velocity of the attachment point at 
x = 0 by Vj and the angular velocity of a body frame fixed m body 1 by u> 


v x = u(o,t): 


- 3 2 U 

“l = '§x - 3t ( °' t)k 


16 


e 


1 


3u 

3x 


(0,t) 


Referring to Figure 2.6, the angle of inclination of the beam with re- 
spect to the x axis is 6 (t) and the vector from the attachment point to 

1 — —1 

the mass center of body 1 is c 1 . If v is the velocity of the mass cen- 
ter of body 1 we can write 


-1 

v_ 


= v 1 + 


x c_ 


where 


c^ = -c^Ccos 0^1 + sin B^j) 


for 


9 1 [ « 1; c 1 s -c 1 (i + 0 1 3) 


Hence 


= [j£ (0 ' t} - c i aHt (0ft) ]^ 

where we have dropped the nonlinear term 00. By straightforward dif- 

li _i 

ferentiation the acceleration of the mass center of body 1 (a ) and the 

. © 

angular acceleration of body 1 (iJ^) are given by 



— ^ (0,t) - c 9 U ? (0,t) j 
3t 3x 3t 


( 2 - 22 ) 


a 3 
3 u 


3x 3t 


(0,t)k 


(2-23) 
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Figure 2.6. Beam and tip body 1 - deformed. 


We now proceed to write dynamic equilibrium equations for body 1. For 
this purpose let S be the force exerted by the beam on the tip body and M 
the moment exerted by the beam on the tip body. 


S 


~2 *3 

3 u . , 3 u 

m — - (0, t) - c, ; 

1 L3t 2 1 3x 3t z 


(0,t) 


(2-24) 


By the 
obtain 


sign conventions of shearing forces in beams S = -El 
the following boundary condition at x = 0 


a 3 
3 u 

3x“ 


so we 


g3 

El (0,t) + m 
3x 


r.2 .3 

3 u , % 3 u 

— — (0,t) - c ; 

,3t 3x 3t" 


(0,t) 


= 0 


(2-25) 


The z component of the time rate of change of the angular momentum of 
body 1 about its mass center is 


J. 9 u ? (0,t) 
3x 3t 


Net moment about mass center 


M - 


x S 
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From beam theory M = El — — (0,t)k and from (2-24) 

3x 2 


■c^ x S = m 


3 2 u 

1 C 1 2 (0 ' 
_3t 


3 3 

t) - c (0,t) k 

3x 3t 


Equating the time rate of change of the angular momentum to the net mo- 
ment we arrive at the second boundary condition at x = 0 


2 2 1 
3 u 3u 12 \ 3u 

El — - ( 0 , t) + m 1 c 1 — — (0,t) - ^m^c 1 + J.J j (O't) 

** ^ 3 x 3 1 


(2-26) 


The boundary conditions at x = Z can be obtained from Eqs. (2-8) and 
(2-9) (recall that m these equations J was the moment of inertia of the 
tip body about the attachment point) 


El — ^ (£,t) - m 
3x 


^ (£,t) + c 5 U - U,t) = 0 

_3t 3x 3t 


(2-27) 


EI 2 (£,t) + m 2 C 2 2 (£,t) + ( J 2 + m 2 C 2 ) 2 


3x 3t 


(2-28) 


We now proceed to the eigenvalue problem associated with the beam bending 
equation and boundary conditions (2-25) - (2-28) . Assuming a solution 
u(x,t) = e la>t <(> (x) we arrive at the eigenvalue problem 


d <f> 4 , . 4 2 . 

— - - a <|> = 0 (a = pw / EI) 

dx 4 
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Boundary conditions at x = 0 


m a 

(0) + — £— [-$(0) + Cl <f>' (0)] = 0 


(2-29) 


m„ c 


♦ " (0) ^ a 4 <}>(0) + y- (0) = 0 (2-30) 


Boundary conditions at x = % 


4 

<J>'"(£) + m 2 y- [<}>(£) + c 2 <j>'U)] = 0 


(2-31) 


m„c„ „ 4 , 

<t>" U) - a 4 <f 1 (£) - y ( 


J 2 + m 2 C 2/^ ' W = 0 


(2-32) 


Orthogonality Relation 

The kinetic energy, T, of the system is given by 


T = j f u 2 (x,t) p dx + 
0 


Im^fO't) - c i (0 ' t} ] 

+ 1 J i[lxTt (0 ' fc) ] + \ m 2 [ li(A ' t) + c 2 3^ < A ' tJ ] 


Expanding u(x,t) m a series of eigenfunctions and assuming the quadratic 

• • 

form in q^ftjq (t) is diagonal we arrive at 
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I 

p / <j> . (x) <j> . (x) dx + m [<(). (0) - c 1 <()! (0) ] [<(> (0) - c_ 4> ! (0) ] 
q i 3 xx xx j xj 

+ (0) 4)^ (0) + m 2 [<j) i (Jl) + c 2 *^(£)][$ 3 W) + (A) ] 

+ J 2 cf>^U)<^(£) = 0 

(2-33) 

where <J>^(x) and (x) are eigenfunctions corresponding to distinct eigen- 
values. 

It can be readily shown that <o = 0 is an eigenvalue of the problem 
corresponding to two linearly independent eigenfunctions: rigid body 

translation, and rigid body rotation. 

Nonzero Bending Modes (a ^ 0) 

The eigenfunction has the form 

<p (x) = A^ sin ax + A 2 cos ax + A^ sinh ax + A^ cosh ax (2-34) 


Applying the boundary conditions at x = 0 (2-29) , (2-30) yields 

^m*c*f) 2 - - m*3A 2 + ^1 + m*c*B 2 jA 2 - m*3A 4 = 0 (2-35) 

(m*c* 2 + J^6 3 A 1 - (l + m £ c £8 2 )A 2 + (m*c* 2 + J*Je 3 A 3 + ^1 - m*c*B 2 ^A 4 = 0 

(2-36) 

where we have introduced the dimensionless tip body parameters 



— and 3 

P* 


al 
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K>* 


Applying the boundary conditions at x = 1 ( 2 - 31 ) , ( 2 - 32 ) yields 


8 sin 8 - ^1 - m*c*8 2 j cos bJa^ + |^m*8 cos 8 + ^1 - m*c*8 2 ^ sm 8^A, 


[• 


+ |m *8 sinh 8 + ^1 + m*c* 8 2 j cosh bJa_ 


+ |m *8 cosh 8 + ^1 + m 2 C 2 ^ 2 ) sin ^ ^] A 4 = 0 

( 2 - 37 ) 


17 2 

j sm 8 + 

/ _ 2 \ 

- [(l + m*c *8 

(J| + »*C* )l 

+ [(j* + m*c * 2 

^ 8 3 sin 8 

- ll + m*c*B 

+ [(l - m*c*B 2 ) 

sinh 8 - { 

[ J 2 + ") c 2 2 ) 6 

+ [i 1 ' m 2 C 2 02 ) 

cosh 8 - | 

h + m 5 ° 2 2 ) B 


i 2 cos bJa^ 
cos 8 ^a 2 
5 cosh bJa^ 


= 0 


( 2 - 38 ) 


The dimensionless tip body parameters m* , c*, J* are defined exactly as 
those for body 1. 

Equations ( 2 - 35 ) - ( 2 - 38 ) constitute a system of 4 homogeneous 
linear equations m A^ , A^, A 3 and A 4> A nontrivial solution exists if 
and only if the coefficient matrix is singular. 

Writing the equations m matrix-vector form with 



[M] A = 0 


Note that 


M 


M(8; m 


1 ' 1 ' 


J l' 


m 


2 ' 2 ' 


J 2> 
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The permissible values of 3 are determined from det[M] = 0 and are func- 
tions of only the six dimensionless tip body parameters even though there 
are nine system parameters. The natural frequencies w are obtained from 

JC 

\ (k = i ' 2 ' •••> 

— (k) — 

Corresponding to 8 , there will be a vector A ^ 0 which in conjunction 
with Eq. (2-34) gives the corresponding eigenfunction. 
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CHAPTER 3 


FORCED VIBRATION OF A CANTILEVERED BEAM WITH TIP BODY 
SUBJECT TO BASE ACCELERATION 


3. 1 Exact Solution 

In Section 2.2 the natural frequencies and mode shapes for a can- 
tilevered beam with tip body were investigated. We now wish to examine 
the more general situation in which the root of the beam (x = 0) is not 
mertially fixed and external forces are present. As noted m Sec- 
tion 2.2 we must be careful that no axial loads (impressed or inertial) 
are acting on the beam if we want to use the simple bending theory. 
Referring to Figure 3.1 let the x-y frame be a body fixed frame attached 
to the left end of the beam (x = 0) and denote by u(x,t) the elastic 
deflection of the beam along the y axis. The beam has a translational 
acceleration a^t) at x = 0 directed along the y axis and the body frame 



Figure 3.1. Cantilevered beam with tip body subject to base acceleration 
and external force and torque. 
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has an inertial angular acceleration ai^ft) perpendicular to the plane of 

motion (along z axis) . Let f (t) denote the external force acting upon 

J? 

the tip body through its center of mass directed along the y axis, and 
g^(t) the external moment on the tip body directed along the z axis. 

The partial differential equation for the elastic displacement 
u(x,t) , is essentially (2-11) modified by D'Alembert's principle. 

a 4 a 2 

El — ^ + p = -pa Q (t) - pxio Q (t) (3-1) 

3x at 

By definition of the body frame we have two geometric boundary conditions 
at x = 0 


u(0,t) = — (0,t) = 0 for all t > 0 (3-2) 

oX — 

The natural boundary conditions at x = i are obtained exactly as 
in Section 2.2 taking into account the effects of a Q (t) and (t) 

a 2 3 2 

El — j (Z, t) + me — ^ (Z,t) 

ax at 

+ j - ■■■ — 2 (Z,t) 

3x at 


= g (t) + cf (t) - mca (t) 

ir ir 


- (mcZ + J)oj (t) 


(3-3) 


a 3 a 2 

EI — | (Z,t) - m — £ (A,t) 

3x at 


- me 3 U r U,t) 
3x at 


-f (t) + ma Q (t) + m (Z + c)a) Q (t) 


(3-4) 
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The boundary value problem (3-1) - (3-4) consists of a nonhomo- 
geneous partial differential equation and time dependent boundary condi- 
tions. A key m solving this problem consists of finding solutions of 
the associated homogeneous problem for which the base acceleration and 

external excitation is zero, i.e. , a. (t) = or (t) — g (t) = f (t) = 0. 

0 0 p p 

This problem has been solved m Section 2.2 where the eigenfunctions 

<}> k (x) and eigenvalues were derived (see Eqs. (2-21) and (2-19)). It 

will prove more convenient to work with the dimensionless eigenfunctions 

S k (n) where n = x/i,. These are essentially given by (2-21) with the 

replacement — -*■ n. The eigenvalue problem for S (rj ) can be obtained from 
X/ K 

(2-12) - (2-15) 


d s k (n) 
dn 4 


- X k S k (Tl) 


0 


s k ( 0 ) 


°; — S k (0) = 0 


dn 


3 S k (Tl) 




“ 

= m*A 

i k 

n=i 

V 1 ’ * <=* It V "> 

n=i 




1 

l_ 

_l 

d 2 , , 

2 S k (n) 

dn 

n=i 

ii 

-‘•V 1 ’ + J *ST s k (I '» 

n=i 


(3-5) 


The eigenvalue problem for S k (n) is m terms of the dimensionless tip 
body parameters m*, J* and c*. The natural frequency of vibration is 
related to the eigenvalue A by 




e^-s 2 

El k 
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The orthogonality condition (2-16) gives the natural inner product for 
the eigenfunctions S^n). We assume henceforth that the eigenfunctions 
are orthonormal with respect to this inner product.* 


1 

/ s (n)s (n) dn + m*s (i)s (i) 
o 1 3 1 3 

+ J*S^(l)S' (1) + m*c*[S i (l)S' (1) + S^(1)S (1)] = 6 


The boundary value problem 
tor notation as 

D[u(x,t)] = 

D 1 [u(x,t) l x=0 = 0; 

D^[u(x,t) ] = f 1 (t); 


(3-1) - (3-4) can be written m opera- 

-pa Q (t) - pxw Q (t) 

[u (x,t) ] x _ 0 = 0 

D 4 [u(x,t)] x=Jl = f 2 (t) 


where the partial differential operators D, D^, D^, D^, are 


D 


El 



+ P 



= 1 ; 


_ 3 _ 

3x 


2 2 3 

3 3 3 

D = El — — + me — j + J J > 

3x 3t 3x 3t 





3x 3t^ 


and 


* 

The details of this normalization are given in Appendix A. 
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f- L (t) = g p (t) + cf (t) - mca Q (t) - (mcd + J)a) Q (t) 

f 2 (t) = -f (t) + ma 0 (t) + m(£ + c)w Q (t) 

We write u(x,t) as the sum 

u(x,t) = v (x / 1 ) + h 1 (x)f 1 (t) + h 2 (x)f (t) 

The approach is to choose the functions h^, h 2 such that the boundary 
conditions on v(x,t) are rendered homogeneous. In order that 
D 1 [v(x,t) ] x=0 = 0 we must have h 1 {0)f 1 (t) + h 2 (0)f 2 (t) = 0 for all t > 0. 
This will be satisfied if 1^(0) = 0 and h 2 (0) = 0. Similarly 

D 2 [v(x,t)] x=Q = 0 if h£(0) = 0 and h£(0) = 0 

Ojlvtx.tn^ - £ l( t) - D 3 [h 1 (x)£ l (t)] x , J! - D 3 [h 2W f 2 (t)J x=t 

D 3 [v ( x,t) ] will be zero if 

[EIh“ (£) - l]f (t) + Elh^ ' {£) f (t) 

+ [mch 1 (£) + Jh^(£)]f (t) 

+ [mch 2 (£) + Jh^(£)]f 2 (t) = 0 (for all t > 0) 

The above condition will be satisfied if 

Elh” (£) - 1 = 0; h 2 ' (£} = 0; 

mch 1 (£) + Jh^(£) = 0; mch 2 (£) + Jh^ (£) = 0 
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Similarly [v (x, t) ] = 0 if the following conditions are met 
' (£) = 0; Elh^ ' ' (£) - 1 = 0; 

h^l) + ch^(£) = 0; h 2 U) + ch£{£) = 0 

In summary, we must find functions h^ (x) and (x) which satisfy 
the respective equations 

h^"(£) =0 h (0) = 0 

h x (0) =0 h^(0) = 0 

h£(0) =0 h^’(£) = 0 

Elh” (£) =1 h 2 (£) = 0 

h 1 (£) = 0 EIh”’(£) = 1 

hj_(£) = 0 h^ (£) = 0 


Clearly, these conditions do not determine h^ and h „ uniquely. For con- 
venience we use fifth degree polynomials. Applying the above conditions 
to a general polynomial of fifth degree we obtain 


2 4 3 

x - j x + 


h 2 (x) 




3 

x 



(3-8) 
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We must still construct the function v(x,t) which m addition to 
satisfying a particular differential equation (considered below) has to 
satisfy the homogeneous boundary conditions 

D 1 (v(x,t) ] x _ 0 = 0 ; D 2 [v(x,t) ] x=0 = 0; 

D 3 tv(x,t)] x=jl = 0; D 4 [v(x,t)] x=£ = 0 

If we formally write 

00 

v(x,t) = Z E q^(t) S^(ri) 


the first two boundary conditions are satisfied for arbitrary q (t) . 

k 

However, since the differential operators and contain time deriva- 
tives as well as spatial derivatives the last two boundary conditions can 
not be met for arbitrary q , indeed constraints amongst these coordinates 
are obtained by demanding that these two boundary conditions are satis- 
fied. Inserting this expansion for v(x,t) into the last two boundary 
conditions (assuming term by term differentiation is valid) and using the 
results (3-5) we obtain the following 

00 C" 2 1 s k' (1) 

J J q k (t) + W t} ] -T~ = 0 (3 “ 9 > 

K — X JC 


- r- 2 i s i" (1 > 

+ a A (t) J x7 ■ 0 

“X JC 


Recall that 


v(x,t) = u(x,t) - h 1 (x)f 1 (t) - h 2 (x)f 2 (t) 
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and 


D[u (x,t) ] = -pa Q (t) - pxa) (t) 


Hence the partial differential equation for v(x,t) is 


9 v 9 v 

El — - + p — — = R (x,t) 0 < x < l (3-11) 

9x 9t 


where 


• IV IV 

R(x,t) = -pa Q (t) - pxu) Q (t) - EIh 1 (x) f 1 (t) - EIh 2 (x) f 2 (t) 

- ph 1 (x)f 1 (t) - ph 2 (x) f 2 (t) (3-12) 

(The superscript "iv" indicates a fourth derivative with respect to the 
argument of the function) . 

Inserting the series expansion for v(x,t) and using the differen- 
tial equation for the eigenfunction S (n) ; (3-11) becomes 

K 

pAz(q k + ^k) S k (n) = R(x ' t} 

g 

We observe that R(0,t) f 0 and — R(0,t) f 0 in general. Therefore the 

o X 

above expansion for R(x,t) cannot hold pomtwise on 0 < x < i. We will 
interpret the above expansion m the sense of convergence m the mean. 
This will suffice for our purposes since we will immediately form an 
integral of both sides with (n) . 


p«, Z (q + 
k=l X * 


o 2 > 
ft, q, 
kk 


1 

/ s (n)s, (n) dn 
o 1 k 


l 

/ R(x,t)S (rj ) dp 
„ i 
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(3-13) 



The orthogonality condition (3-6) can be rewritten with the aid of (3-5) 


as 


1 

/ s (n) s (n) dn 
o 1 * 


(l) 


(l) 


Ik 


- S’ (1) 

1 


s ( 1 ) 

1 


This result m conjunction with (3-13) and observance of the identities 
(3-9), (3-10) uncouples the modal coordinate equations for v(x,t). 


P^(q i + 


1 

/ R(x,t)s (n) dn 

o 1 


(3-14) 


This modal coordinate equation can be expanded as 


q + ft q 

n i 11 


F i (t) + t li f l (t> + t 2i f 2 (t) (l = 1 ' 2 ' 


where 


-F (t) 

i 


= j f s (n) dn • a (t) +/ ns i (n) dn w Q (t) 

* 0 1 0 


+ ^7 / ft v (x) s (n) dn f (t) 

P £ o 1 1 1 


+ ^ / h^ v (x) s (n) dn f 9 (t) 

P * g 2 1 2 


= - -r / h (x) S (n) dn 

H * 0 1 


: 0 = - y f h_(x)s (n) dn 

2i SL 2 i 
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It is possible to make a transformation to a new set of coordinates p^(t) 
such that the differential equations on p^(t) do not contain f (t) and 
f 2 (t). Indeed if we set 

q 1 (t) = P i (t) + t^lt) + t 2i f 2 (t) (3-15) 

The modal coordinate equation (3-14) transforms to 

P i + n 2 J P 1 = P i (t) - t li ^f 1 (t) - t 2i ^f 2 (t) (3-16) 

In terms of the new modal coordinates p^(t) the expression for v(x,t) is 


v(x,t) = 11 P k (t) s k ^^ + A 2 t lk S k( n ^ f l (t ) + 1 E t 2k S k^ T1 ^ f 2 (t) 

k=l k=l k=l 


The beam deflection u(x,t) assumes the form 


u(x,t) = 11 p (t)S, (n) + h- (x) + £ E 

k=l k k 1 k=: 


^h x (x) + i E t lk S k^ n ^ f l 


(t) 


h 2 (x) + a E t 2k S k^ n ^ 
k=l 


f 2 (t) 


with the substitution x = Eqs. (3-7) , (3-8) become 


h 1 (x) 


h 2 (x) 


r 

EI 


EI 


/ 3 2 

U" - 


.3 7 4 5 

4n + — n - n 


) H g i 


(n) 


(3-17) 


/ 1 2 13 14 1 5 \ _ 

[- 6 n + 2 11 ' I 11 + 6 n ) = g 2 


(n) 
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For future reference we also have the expressions 


h^ V (x) 


(84 - 120n) 


,1V. X 

h 2 (x) 


(-12 + 20n) 


Consider the expansion of h^ (x) (g^(n)) in terms of the eigenfunctions 


s, (n) • 

k 


g^n) = c i s i^ n ^ + c 2 S 2^ n ^ + 


taking inner products with S (n) yields 


hence 


c k = <s k (n), g x (n)^ = / s k (n)g 1 (n) dn 


5l (l) = g^(l) = 0 


V x) - - “ fc ik s k (T1) 


Similarly it can be shown that 


h 2 (x) = - z £t 2 k s k (n) 

k=l 
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Hence the beam deflection is expressible as a series m the mode shapes 

s k (n) 


u(x,t) 


£ z p (t)s (n) 

k=l k k 


(3-19) 


Returning to the modal coordinate Eq. (3-16) we have 


p + n p = - -T- / s (n) dn • a (t) - / ns (n) dn m (t) 

ill *■ 0 1 u o 


7 / R* h-(x) - fi 2 h(x)ls (n) dn f (t) 

£ 0 L p 1 i i J i i 


(3-20) 


Now 


v [f h r w - a " f 2 (t > 


i 

/ 

0 


I"— h^ v (x) - fi 2 h (x)]s (n) dn = — — 7 / (21 - 30 n)s (n) dn 

Lp 1 1 1 J 1 p£ 2 0 


X i ) (3 2 .3 7 4 

- “2 MI n - 4n + 2 n 

p£ 0 


n 5 ) S i(n) dn 


(3-21) 
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and 


1 r ? i 

f j^~~ h^ v (x) " i h 2 (x)Js i (n) dn = 


— / (-12 + 20 n)s (n) dn 
p *o 


X 1 


- — I - 

oaM 6 


2 13 

n +-n 


14 1 5\ , , 

2 n + g n te(n) dn 


(3-22) 


We must evaluate four weighted integrals of the eigenfunction. This can 
be readily accomplished through integration by parts and use of the geo- 
metric boundary conditions . 


/ n s (n) dn 
o 1 


f s (" (1) -f s ; ,<i > + f 

iii 


/ n s (n) dn 
o 1 


/ n s (n) dn 
o 1 


-[S'"U) 

1 


s^' ' (1) 


3S“ (1) + 6S^(1) - 6S i (l)] 


(3-23) 


4S^' (1) + 12S^(1) 


24S i (l) + 24 / s i (n) dn 


/ n s (n) dn 
o 1 


x 


i 

s^"(i) - 5 S“ (i) + 2os^(i) - 60s^ (i) + 120/ ns (n) dn 

0 J 


Using these results m (3-21) and (3-22) we obtain 


1 r i 

f ^ h ^ V(x) - fi^h (x) s (n) dn = 

0 LM J pz 


2 S l (1) 


/ h^ v (x) - n \ 2 (x)Js i (n) dn = 


^I s i (1) 


36 



Inserting the expressions for (t) and f (t) and using the above re- 
sults, the Eq. (3-20) can be written m its final form 


p. + Q p 
i 11 


= -u 


a 0 (t) 


3i 


U 4i W 0 (t) + U li 


g (t) 
p 

p£ 3 


+ u 


f p (t) 


2i 


P* 


(3-24) 


where the dimensionless modal parameters are defined by 


u 


lk 


s-(D 


u 2k = S k (D + c*S> (!) 

(3-25) 

1 

u 3k = / S k (n) dn + m*S k (1) + m*c*S^(l) 


1 

u., = / ns (n) dn + m*(l + c*)S (1) + (m*c* + J*)S' (1) 

4JC q K K K 


3.2 Modal Parameter Identities 

The modal parameter identities derived m this section are ob- 
tained by many formal operations and all series expansions are to be 
interpreted in the sense of convergence m the mean. Expanding 
1 = c^S^n) + C 2 S 2 ^ + ^d forming the inner product of both sides 

with s k (n) we fmd 

c k = (s k , 1 y = f S k (n) dn + m*S k (l) + m*c*S k (l) = u 3k 
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Hence 


1 


co 


E u 
k=l 


3 A 1 "’ 


Proceeding m the same fashion we have 


n 


k=l 


U 4k S k (Tl) 


Keeping the orthogonality condition (3-6) m mind we can write 


1 

1 = Z Z u_, u_ / S (n)S (ri) dri 

3i 3] n l j 

1-1 A •* 


also 


m* = Z Z u 3l u 3;] m*S i (l)S =) (1) 


Differentiating the series expansion of 1: 0 = Z u S'(l), 

, OJC K 


0 = Z Z u 3i U 3: J*S^(l)S^ (1) 


Multiplying the series expansion of 1 with the above for 0 



therefore 


0 = E E u,.u 0 ,m*c*|s i (l)S' (1) + S^(1)S (l)j 


i D 


3i 3] 


Adding the above results we obtain 


1 + m* = E u 
k=l 


3k 


(3-26) 


To get an identity on we proceed as above. Firstly 


n = e e u u s (n)s (n) 
4i 4] l 1 

13 


therefore 


also 


- = e e u u / s (n)s (n) dn 
13 J 0 J 


m* = E E u. .u. m*S ( 1 ) S (1) 
43 4i 1 3 

1 3 


Differentiating the expansion for n 


1 = 


k=l 


U 4k S k (1) 


hence 


J* = E E u. u. J*S ! (l)S' (1) 
4i 43 13 
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We also have 


1 - J'vV 11 

k 


Multiplying by the expansion above 


1 = 2 * u 4i u 4 D s i (1)s : (1) 


Noting the symmetry of this expression in 1 , j 


2m*c* = 


z i u 4l u 4;] m * c *( s i (1)s j (1) + S^DS^l)) 


It follows by addition of the above results that 


— + m* + J* + 2m*c* 


We now proceed to get an identity on u^ • u^. From the first expan- 


n = n u 3i u 4 _ ) s i(n) (n) 


Integrating over (0,1) 


- = z z u 3 u 4 / s (n)s (n) dn 

i] 1 3 0 1 : 


S T. u u„ m*S I S 1 
. 3i 4] i J 
i 3 


0 = Z u. S’ (1) 

3i i 

.1 
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and 


1 


Z u. S' (1) 
4] ] 


Multiplying these series 

0 = Z I u, u. S' (l)S' (1) 

3i 4] l j 


1 


Z Z u 
i J 


„ u, 

3i 4j 


S 

i 


(l)S^(l) 


and 


0 


Z 

l 


Z u u. S' (1) S (1) 
3i 4j 1 J 


Hence 


m*c* = Z Z 
i 3 


u. u, m*c* 
3i 4} 


(s^DS'd) + S^DS^Dj 


Adding the above results we obtain 


j + m* + m*c* (3-28) 

We can obtain a simple identity amongst the parameters {u ,} by consid- 

2 

ering the expansion of n over (0,1) in terms of the eigenfunctions, 
writing 

2 

n = c s (ti) + c 2 s 2 (n) + ... 


Z u^, u„, 

, , 3k 4k 
k=l 
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we have 


c k = <n 2 , s k (n)> 


c 


k 


/ n 2 s R (n) dn + m*S k (l) + 2J*S R (1) + m*c*(s^(l) + 2 S r (1)) 

O '" 


Evaluating the integral by the first of relations (3-23) and using the 
boundary conditions (3-5) we obtain 


f u ik s k <r,) 

k=l k 


Differentiating this result and evaluating at n = 1 we obtain 


Z 

k=l 



1 


(3-29) 


Now consider the expansion n 


C 1 S 1 (T1) + c 2 S 2 (n) + • • • 


c 


k 


f n 3 S k (ri) dn + m*S k (l) + 3J*S R (1) + m*c* [S R (1) + 3S R (1) 3 
0 


Evaluating the integral by the second relation in (3-23) and invoking 
(3-5) we obtain 


n 


3 


k=l k 


[s£(l) - s k d)]s k (n) 
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or equivalently 


n 3 = 6 E -j— [Cl + c*)u lk - u 2 k ]s k (n) 

k=l k 


Differentiating this result and evaluating at n = 1 we obtain 


1 . ... „ U lk „ U 2k U lk 

— = (1 + c*)E — E — 

k X k k X k 


The first series is given by (3-29) . Hence 


" U lk U 2k 
k=l X k 


1 

I + c* 


Evaluating the expansion for n and its derivative at n=l we can 


1 1 * 

6 2 C 


l f (1 + c * )u ik s k {1) + l f <c * + =* 2| ' 1 ik s k (1 > 

k k k k 


-5 xr [u 2k (s k (1) + c * s k (1))] 

k k 


1 1 
6 + 2 C 


^ r U lk U 2k „ U 2k 
(1 + C*) E ; E 


k X k 


k X k 


Using (3-30) we obtain 


00 u 


2k 


k=l X k 


1 . _2 
= — + c* + c* 


(3-30) 


write 


(3-31) 
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The model parameter identities derived above will prove extremely useful 
for numerical validation of digital simulations. 


3. 3 Approximate Solution 

In the previous section the dynamics of the excited beam with tip 
body was solved in an exact fashion. The partial differential equation 
of motion was derived along with the time dependent boundary conditions. 

A great deal of labor was expended in constructing a solution which sat- 
isfied the governing equation and all boundary conditions. Indeed, it 
will be recalled, that satisfaction of the two natural boundary condi- 
tions at x = i proved most difficult. The end result, (3-19) , was that 
the structural deformation could be expressed as a series m the mode 
shapes S k (n) - mode shapes for a clamped beam with tip body. Although we 
were able to solve the boundary value problem (3-5), (3-6) for ( n ) ; 
these functions depend upon the parameters m*, c* and J* in a complicated 
fashion. The question arises whether we can expand the deformation m a 
series of simpler functions. Specifically, the mode shapes for a clamped 
beam without tip body suggest themselves, since they will serve the pur- 
pose for all beams with any tip body. Since these later eigenfunctions 
satisfy simple homogeneous boundary conditions at x = 2, , the boundary 
conditions (3-3) , (3-4) will not be satisfied when the structural defor- 
mation is expanded in terms of these modes. This will lead to poor con- 
vergence of the series solution compared with the series solution based 
upon the former eigenfunctions. 

We use the symbol () 0 to indicate a function or parameter based 
upon the simple clamped-free eigenfunctions as opposed to the clamped-tip 
body eigenfunctions. Hence 



X°S° 
k k 


0 


s ° (0) = s° ' (0) = 0 
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0 at n = 1 


(3-32) 


— s° 
an 2 * 


^ 3? - 


dn 


3 k 


These functions enjoy the simple orthonormality condition 


/ s° (n) s (n) dn = 6 
o 1 3 13 


Analogous to (3-25) we have the tip body less modal parameters 


= S k’ (1) 


U 2k * S k (1) 


u 


3k 


= f s°(n) dn 


u 


4k 


f ns°(n) dn 
o k 


(3-33) 


Expanding the structural deformation m terms of these modes 

u(x,t) = a l S°(n)P k (t) (3-34) 

It should be understood that this is only an approximation to 

u(x,t) and that we sum over a finite number of terms. The generalized 

coordinates P k (t) will be determined via Lagrange's Equations for which 

expressions are required for the kinetic and potential energies in terms 

of the p, and p, . 

K K 
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The potential energy in the beam is given by 


V 



dx 


Using the expansion (3-34) , integrating by parts, and invoking (3-32) we 
arrive at 


V(t) 


El , o 2 

2l l Vk (t) 

k 


(3-35) 


If we neglect terms of the order structural deformation x angular rate 
then the kinetic energy of the beam can be written 


1 *2 
T = — p / [V + xu + u(x,t)] dx 

r ^ Q u 


where V (t) is the inertial velocity at x = 0, directed along the y-axis. 
Inserting the expansion (3-34) we arrive at 


T r (t) 


\ pi 3 u)g(t) + i p£ 3 E p 3 (t) + pi 2 2 (V Q (t)u° k + £u 0 (t)u» k )p k (t) 
° k k 


1 .,,2 1 . 2 

+ T p^V Q + T p£ Vo 


(3-36) 


Similarly the kinetic energy of the tip body, T , can be expressed as 


T = ~ m[V + a) Z + u(£,t) ] 2 + \ J 

p 2 0 0 2. 


+ me [Vq + aj Q £ 


[“o + TTS 
+ u(a,t)]^U) Q + 3t g x ( , t ) J 


46 



or 


T p = \ m(v 0 + V )2 + 1 J % + mct V v o + V> 

+ £m(V Q + co q £ + C 0 J Q ) Z U 2k P k {t ^ + [Jt °0 + mc (V 0 + ] Z U lk P k (t) 

k k 

+ br m £ 2 I E u! u’ p,p + ijZEu°u°pp + mc£ E E u! u! p p 
2 2 i 23 * 1*3 2 li 13 * 1*3 2 i 13 * 1*3 

13 ID 13 

(3-37) 

The last quadratic form m (3-37) can be rewritten in symmetric form as 


\ mc£ E Z (uj uj + u° ujjp p 

13 


Adding the expressions (3-36) and (3-37) we obtain the total kinetic 
energy. The virtual work 6 w performed by f and g^ is given by 

6w = f <Su(£,t) + (g + cf ) 6 (u (£,t)) 
p p p x 

Using (3-34) we can write 

6 w = E P 6 p 
1 1 


where 


P 

1 


( Si + 


cu“ ) f + 

li p 


u t g 
11 p 


(3-38) 


Lagranges equations are given by 


d 3 t 3T 

dt 3p k ' 3p k 



(k = 1, 2, 3, ...) 
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Inserting the appropriate expressions for T, V and we arrive at the 
system of differential equations 


•• ET 

[M]p + — [K] 
P* 




where the symmetric positive definite matrices [M] and [K] are given by 


M = m*u® u° + J*u° u° + m*c*(u° u° + u° u? ) + 5 
13 2i 2 j li lj 2i lj 2 j li 




(1, 3 =1, 2 , ) 


[K] = diag(A°, A®, A®, ...) 


and the column vectors , V , and are defined as 


V lk U lk 


(3-40) 


V 


2k 


U 2k + C * U Jk 


V 


3k 


U 3k + m * c * u ik + m * u 2k 


V 


4k 


u®^ + m*(l + c*) U 2k + ( m * c * + 0*) u!^k 


(3-41) 


Associated with the system (3-39) we have the eigenvalue problem 


[K] e = u[M]e 


(3-42) 
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Here the eigenvalues p , p , ... are all positive. Denote by [E] the 

matrix whose columns are the eigenvectors e , e , ... and normalize the 

__ ±. £ 

T 

eigenvectors according to [E ] [M] [E] = [I] . Hence 

[E T ] [K] [E] = diag(p x , P 2 , . . .) 

The equations (3-39) can now be uncoupled by applying the linear trans- 
formation 


P = [E]q 


• El T 

£ + — j diag(p , p , ...)q = [E ]v 


pi 


i&) + iE Mi) 

- 1eT ^3 (t) - IETl ^4“0 


(3-43) 


In Section 3.1 we expanded the structural deformation in terms of the 
eigenfunctions S (H) and arrived directly at the uncoupled system (3-24). 

JC 

Here we have expanded the structural deformation in terms of the eigen- 
functions S°(p) and, as would be expected, the modal coordinate equa- 

JC 

tions (3-39) do not uncouple. A transformation was required to arrive at 
the uncoupled set (3-43) . It is reasonable to expect that as the number 
of retained modes S ° (rj) is increased the finite system (3-43) should 

K 

approach the infinite system (3-24) . More specifically 





J*1 

5 

[E T ]Vi 

-2 

s 

T 

[E T ]V 2 

—3 

s 

[e t ]v 3 

—4 

s 

T 

tE T ]V 4 
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CHAPTER 4 


A MODEL OF THE SHUTTLE ORB ITER WITH DEPLOYED PAYLOAD 


4.1 Motion Equations m Terms of Clamped-Tip Body Eigenfunctions 

In this section we derive equations of motion for the vehicle 
presented in Figure 4.1. A large class of shuttle deployed payloads can 
be approximated in this fashion, specifically, as long slender beams with 
attached tip bodies. The orbiter and tip body are assumed rigid while 
the beam is allowed to undergo small, elastic transverse bending. The 
attachment point of the beam to the orbiter is located arbitrarily with 
respect to the mass center of the orbiter, but in keeping with other 
sections of this report, the mass center of the tip body is located along 
the tip tangent of the beam. (It should be noted that the difficulty 
discussed in Section 2.2 is still present due to axial loads exerted by 
the orbiter on the beam. This complicating affect will be ignored m the 
present analysis.) In order to fully understand the interaction between 
the flight control system and the flexible body dynamics, all rigid body 
motion and bending will be restricted to the orbiter pitch plane. Hence 
there are essentially three rigid body degrees of freedom — two transla- 
tional and one rotational, and an infinite number of elastic degrees of 
freedom. Since the flight control system is primarily concerned with 
orbiter attitude control we analytically eliminate the translational 
coordinates so that the final set of differential equations of motion 
only involve orbiter attitude and elastic degrees of freedom. 

The attitude control system exerts a net force F„ and moment on 

1 - 0-0 

the orbiter. There is an external torque g acting on the tip body and 

-p 

external force f applied at its mass center perpendicular to the neutral 
“P 

axis of the beam. 
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ORBITER BODY FRAME FIXED AT ORBITER MASS CENTER (® Q ) 

BODY FRAME AT 0 CANTILEVERED AT BEAM TIP, PARALLEL TO 
' Xj AXIS PASSES THROUGH TIP BODY MASS CENTER <© ) 

e PITCH ANGLE OF ORBITER MEASURED FROM Xj TO x Q 

Figure 4.1. Planar dynamics model of Shuttle Orbiter 
with deployed payload. 


Basic Mass and Geometric Parameters 
m Q - orbiter mass 

I - orbiter moment of inertia about its mass center perpendicular to 
plane of motion 

a - vector from orbiter mass center to beam attachment point 
p - mass per unit length of beam 
Z - beam length 
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c - tip body mass center offset 


m - tip body mass 

I t “ tip body moment of inertia about its mass center perpendicular to 
plane of motion 


It will prove convenient to define the following quantities, 
m. = p£ + m 


m = hIq + m^ 


m 0 /m 


m /m 


pi /2 + m (SL + c) 

m. 


mass of beam plus tip body 
total vehicle mass 

dimensionless mass ratios 


(4-1) 


location of mass center of beam + tip body 
relative to 0 under no deformation 

(4-2) 


Primordial Motion Equations 

Figure 4.2 is the free body diagram of two sub-bodies associated 
with the vehicle being analyzed: the orbiter, and the composite body of 

beam + tip body. F and G are respectively the interbody force and moment 
exerted by the beam on the orbiter at the attachment point 0. 

We have the translational equations 


m„a° 

0 -® 


Iq + I 


1 

m„ 

1 -© 


f - F 

-p 


(4-3) 


where a° and a are the inertial accelerations of the orbiter mass center 
— © -© 

and beam + tip body mass center respectively. The corresponding rota- 
tional equations are 
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Figure 4.2. Free body diagrams. 


3 ,0 
dt 


G + G + a, ^ F 
— 0 — — - 


(4-4) 


^h 1 
dt — © 


-G+g + bxF+(£+c-b+6)xf 
— -p — — — — — — p — p 


(4-5) 

where h°, h?" are the angular momenta of the orbiter and beam + tip body 
■ G ■ © 

about their respective mass centers. The vector b locates the mass cen- 
ter of the deformed beam + tip body relative to 0 and is illustrated 
along with 6 and c in Figure 4-3. 


If denotes the inertial acceleration of the attachment point 0 
then we have the simple conditions 




Figure 4-3. Deformed beam + tip body - vector geometry. 


-^0 


_o , d 

-• at 2 - 


i 


+ 


dt 


2 - 


(4-6) 


Equation (4-5) can be rewritten m terms of h^ - the angular momentum of 
the beam + tip body about 0. Using the second of results (4-3) and (4-6) 
and the translation theorem for angular momentum we arrive at the more 
useful form 


A. w 1 

dt —0 


-G - 


m^b x 


a„ + g + 

— 0 — p 


(£ + c + 


fi ) 

“P 


f 

“P 


(4-7) 


Ignoring completely the affect of axial acceleration on the beam + 
tip body, the elastic deformation of the system in Figure 4.3 is identi- 
cal to the problem considered m Section 3.1. Indeed, if u(x,t) denotes 
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the elastic deflection of the beam along the y axis of the governing 
partial differential equation is given by (3-1) where a Q (t) is now the 
component of aQ along the y axis and u) Q (t) is the inertial angular accel- 
eration of perpendicular to the plane of motion. By definition of the 
body frame & , the geometric boundary conditions (3-2) apply. The natu- 
ral boundary conditions are given by (3-3) and (3-4) with the following 
replacements 


m -*■ m 

2 

J -*■ It + m^c 


It was shown that u(x,t) = % Z p (t) S, (B) with the modal coordinate 

k=l K 

equations given by (3-24) . 

Elimination of Interbody Force F from Rotational Equations 

Adding the translational equations (4-3) and invoking (4-6) we can 

show 


^0 


= m, 


dt 


2 £ - ^1 


d , -1, 

— £ + m + V 


dt 


(4-8) 


Inserting this expression into the first of (4-3) we obtain 


F 


-my 0 y 


>( 


a + 


dt 



^ + 


y O^P 


(4-9) 


which expresses the interbody force m terms of kinematic quantities and 
the external forces on the vehicle. 
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The orbiter rotational motion equation (4-4) can now be written 


— h° 

dt Ab 


Gq + G - n\y Q M 1 a x 



(a + b) 


y x a x Fq + 


u a x 


f 

-P 


(4-10) 


and the rotational motion equation (4-7) becomes 


jL >1 

dt — 0 


d 

-G - mi^b x — - (y^a - y^b) - y^ x Fq + g p 
dt 

+ (£ + c + 6 - y.,b) x f 

- - -P 1- -P 

(4-11) 


Angular Momentum Calculations 

Let i, : , k be unit vectors along the x, y and z axes of ^ (^) 
and let the plane of motion be the xy plane. If 9 is the orbiter pitch 
rate then 


= I Q 0k + ( )i + ( )J (4-12) 

If u(x,t) is the elastic deflection of the beam and we neglect the term 
9u then the angular momentum of the beam about 0 has a nonzero component 
along z given by 




/ px It- (x,t) dx + p£ 3 9 

0 3t 3 


The z component of the angular momentum of the tip body about 0 is 


c t[® + arts <*■«] + “t[‘i 


+ 5 + c) x — (2, + 5 + 

— p — dt — — p 


4 
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We have 


£ + <$ + c 

- -IP - 


(c + l)i + [u(£,t) + c (£,t)J' 


From this point onward terms of the order structural deflection x angular 
rate will be dropped as well as any derivatives of such quantities. It 
then can be shown that the z component of the angular momentum of the tip 
body about 0 is given by 


\ 2 1 r . 3 3 2 

C. 9 + r~ y (£,t) + ra (c + £) (c + £) 0 + ^ (£,t) + c ~r ~ (£,t) 

t[ 3t 3x J t |_ dt dt dX 


Adding these two contributions it follows that 


4 - i ■ V + / px I? lx,t) dx + r t wh a '-> 
- V c + 11 [li 4 <= irk (t - t) ] 


(4-13) 


where is the inertia of the beam + tip body about 0 m the undeformed 
state 


13 2 

J Q = — p£ + I t + m t (c + £) 


(4-14) 


Differentiating (4-12) , (4-13) it follows that 


(iU) 


k . i o e 


(4-15) 
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Z 2 .3 

k = J 0 + / px — ^ (x,t) dx + I — — ^ — U,t) 
0 0 3t 3t 3x 


+ m (£, + c) 


U,t) + c — 3 „— a,t) 


at 


at 2 ax 


(4-16) 


Vector Geometry 

The motion equations (4-10) , (4-11) call for various geometric 
vectors and their time derivatives. Recall that b is the position vector 
of the mass center of beam + tip body relative to 0. 


b = 


b l + — 
1 m. 


u (x , t ) p dx + m^u ( & , t ) + c (&,t)J 


Using the series expansion for u(x,t) and definition (3-25) for the modal 
parameter u 3k we can show 


b = b i + u (t) d 
— 1 c 


(4-17) 


where u (t) , the shift in the mass center of the beam + tip body due to 
c 

elastic deformation, is given by 


u c (t) 



E u 
k=l 


3k 


P k (t) 


It follows that 


_d_ 

dt 


b 


-u (t) 0i + (b. 0 + u ) j 
c 1 c 


s (b 0 + u ) j 
l c 


(since u is of order structural deflection) 
c 


(4-18) 
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( 4 - 19 ) 


d b = -b 0 2 x + (b 6 + u c ) j 


dt 


2 - 


a = 



( 4 - 20 ) 


It follows that 



( 4 - 21 ) 


a x 



a x 



(a 2 + a 2 )ek 

( a i^c + a l b l® + a 2 b l® 2 )^ 


b x 



i b lV - b l a 2® 



b x 



(Vc + b l’ 9 ’)' 


( 4 - 22 ) 
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Writing EV 
— 0 


m 


a x 




a x f 
- -p 


b 




(£ + c + 6 

— - -p 



= (a 1 F Q2 (t) - 


= a l f p (t)k 



a„F„„ 
2 01 


(t) )k 


= ( bl F 02 (t) - U c (t)F 01 (t))k 

y b) x f = a + C - M _b.)f (t)k 
1— — p lip 


(4-23) 


Rotational Motion Equations - Expanded Form 

Using the above results with (4-15) allows us to write (4-10) m 
the form 

[l Q + my^a \ 

+ y i a 2 F 01 (t) + y 0 a i f p (t) " mtJ 0 P l a 2 b l® 2 

(4-24) 

Here G Q (t) is the component of perpendicular to the plane of motion 
(only nonzero component by assumption) and (t) is the corresponding 
component of G. 


+ a 2 + a 


l b l)_ 


+ mM 0 y i a i U c = G 0 (t) + G z (t) ' y i a i F 02 (t) 
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Similarly with the aid of (4-16) we can write (4-11) in the form 


(j Q + ^ 0 U 1 b 1 a 1 - my^b^Q 


1 a 2 a 3 

+ / px — j (x,t) dx + I — r-^ 

0 3t 3t 3x 


a,t) 


+ m ^ ( H + c ) 


Mr U,t> + C ■ ■ \ u - (£ ,t) 


3t 


3t 2 3x 


- my b u 
lie 


-G (t) - y b F (t) + y u (t) F 
Z 1 1 02 1 C 

+ g p (t) + (H + c - y 1 b 1 )f (t) 


01 


(t) 


+ my 0 y 1 b l a 2 e 


(4-25) 


Using the series expansion for u(x,t) and definition (3-25) for the modal 
parameter u^ k we can show that 


/ px Mr (x,t) dx + I — — (£,t) 
0 3t 3t 3x 


+ m^_(J?, + c) 


3 2 u 

3t 2 


(£,t) + c 


a 3 
3 u 

3t 2 3x 


U,t) 


= p * = Wk (t) 


Introducing this result into (4-25) and adding with (4-24) eliminates the 
interbody torque G ^ (t) 
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0 


+ J Q + my. 


ij ,j i ( a i + 4 + 2a i b i) - ”> , i b i] 5 


+ pt £ (•‘ 0 a 1 » 3k + 2 >V - 

k=l 


G 0 (t> - V*! + b l )P 02 (t) 


+ V a 2 + u o it))F 01 (t) + 3 p lt) 


+ (pa + l + c - y b )f (t) 
0 1 lip 

(4-26) 


As discussed above, the modal coordinates p (t) are governed by 

K 

• •• 

(3-24) with u)q (t) = 0 and a^(t) being the component of a^ along the 
y axis. From (4-8) with (4-19) and (4-21) 


no .. 1 -2 

a 0 (t > - <Vl - W 6 - — l Vl + » (P 02 + V ' “o*2 e 

k 


The modal coordinate differential equation assumes the form 


("o -r u 3l ■ "l ~t u 3i + U 4 J S 


+ Z ( 6 

k=l 


(<S - — u u 

\ ik m 3i 3k/ 


g (t) f (t) 

u — + u 

li _„3 2i .,2 


pi' 


pl 


u a 

~ (f + f ) + u, y n ~ 0 2 - fi 2 p (t) 
m£ 02 p 3i 0 l ii 


(i = 1, 2, 3, ...) 


(4-27) 


Note the coupling of the modal coordinate equation with all the elastic 
coordinates as well as orbiter pitch. Introduce the following notation 
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I o + J o , , ( a i a 2 ^ a i b l\ , , 2,2,„; 

~3 + {m/pl) Vi ( 72 + 71 + 2 T * X j - (m/p *> 

p X/ ' X/ X/ * 


2 2 


00 


a b 

a ik - P 0 T u 3k + U 4k " "l T u 3k 1.2,3, ...) 


Equations (4-26) , (4-27) can now be written 


(4-28) 


a oo 8 + “Vk 

k=l 


G 0 (t) 

P* 3 


'a, + b,\P 02 (t) 


- ^ (^) 


Pi 


{a 2 + V F 01 (t) g p (t) 

+ y, • — + — 


1 1 


p£ 2 p£ 3 


('»T 


b,\ f (t) 


+ 1 + c* - y 


11/ 


Pi 


(4-29) 


a. 0 + £ (5 

11 k-P 115 


(3 - — u u ) 

\ lk m 3i 3k/ 


g (t) f (t) 

P k = U li “TT + U 2i T7T 


p£' 


pr 


if (F 02 (t) + f P (t)) 


+ U 3i P 0 f ^ " n i P i (t) 


( 1 — I, 2/ 3/ « • • ) 


(4-30) 


Equations (4-29) and (4-30) are the final set of differential equations 

of motion involving orbiter pitch and elastic coordinates of the beam. 

The coefficients of the accelerations are constant but (4-29) contains a 

time dependent coefficient through the term F (t)u (t) while (4-30) has 

U J- c 

a nonlinear term in 0 - 
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Matrix-Vector Form of Motion Equations 
Define 


*o (t > 


0 (t> 1‘ 1 + M 

G + c * + “o T ‘ “1 t) 


F 0 2 (t) . *2 F 01 (t) . %«=> 
„2 + U 1 l „2 + . 3 


p£ 


PA' 


f (t) 
_P 

pjl 2 


^ (t) 
1 


= u 


ll 


9 p (t) 


pt 


£ <t> u 

F + "21 ^--ST (P 02 (t » + f p (t ” 

pH 


(i = 1, 2, 3, ...) 


The generalized forces ^(t) and are functions only of the 

forces and moments acting upon the vehicle. Let 


x = 


- ( 9 / P 2 ' P 3 ' • • • ^ > & ~ Q^q ' & p ' ^2 ' ^ 3 ' * • ■ ^ 


N = 


u c (t) F 01 (tf 


0 

M 1 £ 2 



PA 


ii 



31 

0 

a 2 



+ u — 

u 


o 

32 

0 






0 

• 

33 

— 


• 


e 2 (t) 


(4-31) 

external 


(4-32) 
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We have the symmetric generalized mass matrix A 


A = 


00 


11 


12 


13 


a ll 


a i2 


a i3 


(l- 

M 

2 \ 
u.,, ) 

_ m 



_ M 

u 

u 

\ 

m 

31/ 

m 

31 

32 

m 

31 

33 

_ £l 

u „ 


(i- 


2 \ 
u „ / 

_ £l 

u 

u 

m 

32 

31 


m 

32/ 

m 

32 

33 

_ £i 

u„„ 

u 

_ 

u 

u 

(l - 


2 

u 

m 

33 

31 

m 

33 

32 

\ 

m 

33 


and stiffness matrix 


[K] 


( 2 2 2 \ 

0, fi , ft 2 , 0 3 , J 


In terms of the above notation Eqs. (4-29) , (4-30) are 


2 - _ _ _ _ 
[A] — ~ + [K] x = &{t) + N (t; x,x) 
dt 


(4-33) 


(4-34) 


4.2 Motion Equations in Terms of Clamped-Free Eigenfunctions 

We treat here the same vehicle as in Section 4.1; the only differ- 
ence being the set of eigenfunctions used to expand the elastic displace- 
ment m the beam. Here we use the tip-bodyless eigenfunctions S°(n) 
described in Section 3.3. For the same reasons as given there we expect 
that, for the same degree of truncation, the results will be inferior to 
those obtained by use of the motion equations of the previous section. 

The results of the present section in conjunction with our previous anal- 
ysis can serve as a basis for numerical confirmation of this statement. 
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Throughout this section it will be understood that the elastic 
deflection m the beam is expanded m terms of the set {S^(rt)} 


u(x,t) = 1 Z s£(n)P k (t) 


No confusion should arise between the modal coordinates p, (t) and those 

k 

appearing in Section 4.1. Virtually all results previously obtained 
remain intact except in those instances where explicit use was made of 
the series expansion for u(x,t) . We list sequentially those equations m 
Section 4.1 which require modification due to a switch of the modal basis 
from {s k (n)} to {s k (n) >. 

(1) Equation (4-18) is replaced by 
P * 2 

u (t) = E v ,.p. (f) (V is given by Eq. (3-41)) 

c m, , 3k K 3k 

1 k 


(2) Equation (4-26) - second term on left hand side is replaced 
by 


pS. 3 Z 


("0 T v 3k + V 4k - “l T V 3kX <l> 


(3) The modal coordinates are now governed by Eq. (3-39) with 


• • A 0 •• 1 *2 

a (t) = (y a — y b ) 0 — ZV p + — (F + f ) — y a 0 

CT ' ' M 0 1 l r m r_ 3k P k m v 02 p ; M 0 2 


and = 0. Thus Equation (4-27) is replaced by 
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4i 


a i b i V 

+ p o T V 3i - P 1 T V 3i/ 

+ E kk - f v 3 l V 3 kk 


f 

v — + V — 2 — (F + f ) 

li 3 2i 2 ml u 02 p ; 
p£ p£ 


+ V 3 i P 0 f ^ - ~4 X k (t) 

P«- 


where M is given by Eq. (3-40) 

IK 

(4) Equation (4-28) is replaced by 


d l b l 

a lk P 0 £ V 3k + V 4k P 1 £ V 3k (k 1 ' 2 ' 3 ' 


(5) Equation (4-30) is replaced by 


p£ 


a 0 + Z M . - ^ V, V„ 
li \ lk m 3i 3k/ k 

X 




g f V, 

V — 2- + v — 2 i. (F + f ) 

li .3 2i 2 m£ v 02 V 

p£ p£ 


+ Vo l X (t > 

p£ 


(with a given by (4) above) . 

(6) The second definition (4-31) is replaced by 


«^(t) 


g f V, 

E v -. + V. — T (F + f ) 

li .,3 2i _ „ 2 m£ 02 p 


pi' 


Pi 


(7) In Eq. (4-32) use the expression for u^(t) given by (1) 
above and replace second vector with 


’oT i!|ll “ V 31 V 32 V 33 V 34 •••> T 
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( 8 ) 


Definition (4-33) is replaced with 


00 


11 


12 


13 


A = 


l ll l M ll 


l 12 \ M 21 


13 


(m,, ~ ~ 

\ 31 m 


_ 2k 

V V ) 

( M 

- 2k 

V V „ 1 

(M 

_ 2k 

V V 

m 

31 31/ 

l 12 

m 

31 32/ 

\ 13 

m 

31 33 

- 2k 

V V ) 

[ M 

2k 

V ) 

(m 

- 2k 


m 

32 31/ 

\ 22 

m 

32 32/ 

\ 23 

m 

32 33 

_ 2k 

V V ) 

(m 

_ 2k 

V V ) 


- 2k 


m 

33 31/ 

\ 32 

m 

33 32/ 

\ 33 

m 

33 33 










(9) Definition (4-34) is replaced by 


[K] = — diag (0 , Xj, X°, X°, ...) 

pi 


Note: these are the squares of the natural fre- 

quencies of a uniform clamped beam without 
tip body. 

4. 3 Differential Equations of Motion - First Order Form 

In the last two sections equations of motion were derived for the 
vehicle depicted in Figure 4.1. These equations were of the form of a 
system of second order differential equations and are not directly suita- 

d 2 x 

ble for digital implementation. Since these equations are linear in — - 

dt 

they can be transformed into an equivalent set of first order differen- 
tial equations. In effecting this transformation the "system natural 
frequencies" reveal themselves and the option of further modal truncation 
(beyond the truncation of elastic degrees of freedom in the beam) ' can be 
adopted. 
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Our motion equations can be written 


[A] - X + ~~ [B] x = t ) + N ( t; x, x) 

dt p l 


(4-35) 


where 


[B] = diag (0, X , X y ... X^) 

If we follow the procedure in Section 4.1 the X are the eigenvalues of 
the problem for a clamped beam with tip body (Equation 3-5) while if we 
follow the procedure in Section 4.2 they are the eigenvalues appropriate 
to a clamped beam without tip body (Equation 3-32) . 

The free vibration of the system associated with (4-35) leads to 
the eigenvalue problem 


[B]V = y [A] V 

s 

and the "system natural frequencies" a are given by 


(4-36) 


s 


0) 



[B] is symmetric, positive-semidefinite and [A] is symmetric, positive 

definite. All eigenvalues are therefore nonnegative. Specifically y = 0 

T 

is a simple eigenvalue with eigenvector (1, 0, 0, ..., 0) corresponding 
to rigid body rotation. This eigenvalue problem has a full set of line- 
arly independent eigenvectors , ..., V^ n+1 ^ and we can always 

arrange that they be orthonormal relative to [A] , i.e. 


V (l)T [A]V (;3) 


6 

13 
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Define 


[Z] - [v a > v ' 21 v l3 » 


V (n+ll I 


If we make the linear transformation x = [Z]q in (4-35) we 


at 


pT T • 

q + —■ diagUi u • • • V* >3. = [Z ] (^(t) + N(t; x, x) ) 

2 — .4 12 n+1 -*• — — — — 

dt p£ 


Let 


y (t) = 


q(t) 

£(t) 


and label the eigenvalues so that = 0. 

The system (4-38) can now be written as 



y n + 2 


Y n + 3 


y n + 4 


y 2n+2 

m t>T 

[Z ] + N) (0, 

p£ 

P 2 y 2' y 3 y 3' Vl y n+1 )T 




(4-37) 

arrive 


(4-38) 


(4-39) 
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CHAPTER 5 


FORTRAN PROGRAMS 


Two separate FORTRAN computer programs have been created for the 
purpose of numerically implementing the motion equations of Sections 4.1 
and 4.3. Complete listings of these programs accompanied by annotated 
sample input and output data are provided m Appendices B and C. Each is 
extensively commented throughout and in most instances the FORTRAN vari- 
able names are mnemonically similar to corresponding analytical quanti- 
ties. Where relevant, reference is made to specific equations of the 
report. 

The program of Appendix B computes the eigenvalues, modal param- 
eters and modal parameter identities of a cantilevered beam with tip 
body. Note that these computations are performed in (IBM) quadruple 
precision due to the numerical sensitivity of the transcendental expres- 
sions involved. The roots of the characteristic equation (Eq. 2-19) are 
estimated by incrementing the parameter 3 and searching for sign changes 
in the left-hand-side. Root estimates are then improved by Newton- 
Raphson iteration. The modal parameters of Eq. 3-25 are computed and 
together with the roots of Eq. 2-19 (raised to the fourth power) written 
to a disc data set to be used by the motion equation program (described 
below) . The partial sums and asymptotic values of the modal parameter 
identities (Eqs. 3-26 thru 3-31) are evaluated and output. This informa- 
tion serves to check the calculations of the modal parameters and can 
indicate "missed" roots of the characteristic equation by showing poor 
convergence . 
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The program of Appendix C numerically synthesizes and integrates 
the motion equations for the planar dynamics of an elastic beam with 
rigid bodies at each end. The formulation implemented is that for which 
the beam deflection is expanded as a series of eigenfunctions appropriate 
to a clamped beam with tip body (as per Section 4.1) . The eigenvalues 
and modal parameters generated by the above program are read and checked 
for consistency with the NAMELIST input data. Note that the initial 
conditions on the attitude angle and its rate may be specified arbitrar- 
ily as input, while those of the modal coordinates and their time deriva- 
tives are internally set to zero. The algebraic eigenvalue problem of 
Eq. 4-36 is solved via the double precision IMSL subroutine EIGZS. The 
"system natural frequencies" are output and are helpful for selection of 
an integration time step. External forces and torques on the "Orbiter" 
and tip body are specified through subroutines ORBFOR and TPBFOR respec- 
tively. The motion equations m first order form, Eq. 4-39 are inte- 
grated using the Adams Method with third order differences. 
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APPENDIX A 


EIGENFUNCTION NORMALIZATION FOR CANTILEVERED BEAM WITH TIP BODY 


In Section 3.1 extensive use was made of the eigenfunctions S (n) 

K 

defined by the boundary value problem (3-5) . In this appendix we formu- 
late a procedure to obtain these normalized eigenfunctions from a set of 
non-normalized eigenfunctions. Certain definite integrals involving the 
eigenfunctions are evaluated to facilitate computation of the modal pa- 
rameters given m Eqs. (3-25) . 

Recall that for the case of a simple beam with no tip body the 
normalized clamped-free eigenfunctions are given by 


sm 3^ - smh 3^ 

V n) = cosh 3 k n - cos B k n + ^rrr^hT (sinh V - sin 6 k n) 


where 


is a root of 


cosh 3 cos 3+1 =0 

and the normalization condition is 

1 2 

f y <n) dn = l 
0 k 

To obtain a set of eigenfunctions for the case of a beam with tip 
body we can, m the discussion preceding Eqs. (2-17) , (2-18) , choose the 
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constant c, =1. (So c„ = -1) . Invoking the boundary condition (2-18) 
4 2 

the resulting unnormalized eigenfunction ip (g) is given by 


^^(ri) = cosh - cos 3 k n + Y^lsmh B k n - sm B k n) (A.l) 


where now B k are the roots of Eq. (2-19) and the parameter y k is given by 


2 3 

m*c*B (cosh B - cos B) + J*B (sm B + smh B) - cos B - cosh 3 
Y ~ 2 3 

m*c*B (sm B - sinh B) + J*B (cos B - cosh 3) + sm B + smh B 

where subscript k is deleted for clarity. (Note that for m* = J* = 0 
this eigenfunction is already normalized) . 

Using the eigenfunction (A.l) we readily evaluate the following 
integrals 


/ i > 1 2 (n) dn 
o 


1 + — (sinh 8 cosh 8 + sm 8 cos 8) 
28 


- — (sinh 8 cos 8 + cosh 8 sin 8) 
8 


Y 

+ — [sinh B cosh 3 - sin 3 cos 
23 


- 2(sin B cosh B ~ cos 3 sinh 3)1 

v 2 

+ — (sm B - smh 8) 

B 

(A. 2) 

1 1 

/ ^ (g) dg = — [smh B - sin B + y(cosh B + cos 8 - 2) ] (A. 3) 

o p 
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1 

/ ruMn) dr > 
o 


\ [smh 8 - sin B + Y (cosh B + cos B) ] 
B 


— [2 - cosh B - cos 
8 


We can now evaluate the square of the norm of 
with respect to the inner product m (3-6) 


P 


k 


1 

/ 

0 


^(n) dn + (i) + j* 



The normalized eigenfunctions are given by 


s. (n) 

k 


-£= V n) 


- y (sin B + smh B) ] 

(A. 4) 

the eigenfunction (n) 

+ 2m*c*iK (1) 4i'(l) 

K K 

(A. 5) 


75 



APPENDIX B 


"EIGENVALUE "-MDDAL PARAMETER 
FORTRAN PROGRAM LISTING 


ISN 0002 
ISN 0003 
ISN 0004 
ISN 0005 


ISN OOOG 
ISN 0007 
ISN 0008 
ISN 0009 
ISN 0011 
ISN 0012 
ISN 0013 
ISN 0014 
ISN 0015 
ISN 0016 
ISN 0017 
ISN 0018 
ISN 0019 
ISN 0020 


ISN 0021 
ISN 0022 
ISN 0023 
ISN 0024 


C ********************************************************************* *00000100 

C * THIS PROGRAM COMPUTES THE EIGENVALUES. MODAL PARAMETERS AND *00000110 
C * MODAL PARAMETER IDENTITIES OF A CLAMPED BEAM WITH TIP BODY *00000120 
C * THE EIGENVALUES AND MODAL PARAMETERS WHICH ARE NECESSARY INPUT *00000130 
C * TO THE PLANAR DYNAMICS PROGRAM ARE WRITTEN TO A DISC FILE *00000140 
C * ( WRITTEN BY JOEL STORCH 8 STEPHEN GATES C S D L BASED UPON *00000200 
C * C S D L REPORT # R-1629 MAY 1983 ) *00000300 
q ********************************************************************* *00000400 


NOTE ARRAYS DIMENSIONED FOR A MAXIMUM OF 50 EIGENVALUES 

IMPLICIT REAL* 16(A-H,0-Z) 

REAL* 16 LAM.MSTAR, JSTAR 

DIMENSION LAM( 50) , U 1 ( 50) ,U2( 50) , U3( 50) , U4( 50) 

NAMELIST /INPUT/ MSTAR . JSTAR , CSTAR , NF 


00000500 
00000510 
00000520 
00000600 
00000700 
00000900 
00000910 
00000918 

********************** input - OUTPUT FILES ************************ *00000926 
FILE #5 NAMELIST INPUT FILE 00000934 

FILE #6 PRINTED OUTPUT FILE 00000942 

FILE #8 DISC OUTPUT FILE FOR EIGENVALUES 6 MODAL PARAMETERS 00000950 

00000958 

**************** DESCRIPTION OF /INPUT/ LIST ITEMS ******* ***********00000966 

00000974 

"MSTAR” MASS RATIO (OF EQ 2-19) 00000982 

"JSTAR" INERTIA RATIO (OF EQ 2-19) 00000983 

"CSTAR" OFFSET RATIO (OF EO 2-19) 00000984 

"NF“ NUMBER OF EIGENVALUES & MODAL PARAMETERS TO BE COMPUTED 00000985 

********************************************************************* 00000986 

00000990 

PROMPT FOR NAMELIST INPUT DATA 00000993 

00000996 
00001000 
00001100 
00001200 
00001300 
00001400 
00001500 
00001600 
00001610 
00001700 
00001800 
00001810 
00001820 
00001830 
00001900 


WRITE(5, 1) 

FORMAT ( 1H ,5X. 'INPUT DATA',/) 

READ(5, INPUT) 

IF(NF LE 50) GO TO 3 

WRITE(6,2) NF 

STOP 

FORMAT ( 1H0, 5X , 13 , ' MODES REOUESTED MAXIMUM ALLOWABLE IS 50') 
WRITE(6,6) 

WRITE(6,4) MSTAR, JSTAR.CSTAR 

FORMAT ( 1H0.5X. 'MSTAR 3 ' ,F8 4, 4X , ' JST AR* ' . F8 4 . 4X . ' CSTAR* ' , F8 4) 

WRITE(6,5) NF 

FORMAT ( 1H ,5X.'NF='.I3) 

F0RMAT(1H1,1SX,'0ATA FROM NAMELIST INPUT') 

CALL EIGEN(NF, MSTAR, JSTAR, CST AR , LAM , U1 , U2 , U3 , U4 ) 


00001910 

OUTPUT EIGENVALUES AND MODAL PARAMETERS (NOTE ARRAY LAM CONTAINS THE00001920 
ROOTS OF EO 2-19 RAISED TO THE FOURTH POWER) 00001930 

00001940 

WRITE(6, 13) 00001950 

WRITE (6, 11) 00002000 

WRITE (8 ) NF, MSTAR, JSTAR.CSTAR 00002100 

DO 10 N= 1 , NF 00002200 
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ISN 

0025 


WRITE(8) LAM(N),U1(N),U2(N),U3(N).U4(N) 

00002300 

ISN 

0026 

10 

WRITE(6, 12) N,LAM(N),U1(N),U2(N),U3(N).U4(N) 

00002400 

ISN 

0027 

11 

FORMAT ( 1HO,T4,'N',T12, 'LAM' ,T27, 'U1 ' ,T40, 'U2' ,T49, 'U3 ' , T66 , 'U4' ) 

00002500 

ISN 

0028 

12 

FORMAT! 1H0.T3, 12. T8.G12 5.T24.F9 4. T36. F9 . 4.T47 , F9 5.T60.G13 5) 

00002700 

ISN 

0029 

13 

FORMAT! 1H0./. 17X. '"EIGENVALUES" & MODAL PARAMETERS') 

00002710 



C 


00002800 



C 

MODAL PARAMETER IDENTITIES 

00002900 



C 


00003000 

ISN 

0030 


WRITE (6 . 20) 

00003 100 

ISN 

0031 

20 

FORMAT! 1H1.35X, 'MOOAL PARAMETER IDENTITIES') 

00003200 

ISN 

0032 


WRITE ( 6,23) 

00003210 

ISN 

0033 


SUM 1=0 000 

00003300 

ISN 

0034 


SUM2'0 000 

00003400 

ISN 

0035 


SUM3=0.0Q0 

00003500 

ISN 

0036 


SUM4=0 000 

00003600 

ISN 

0037 


SUMS'O OOO 

00003700 

ISN 

0038 


SUM6=0 OOO 

00003800 

ISN 

0039 


00 22 N= 1 , NF 

00003900 

ISN 

0040 


SUM1=SUM1+U3(N)**2 

00004000 

ISN 

0041 


SUM2 = SUM2HI4(N)**2 

00004100 

ISN 

0042 


SUM3=SUM3+U3(N)*U4(N) 

00004200 

ISN 

0043 


SUM4'SUM4+U1(N)**2/LAM(N) 

00004300 

ISN 

0044 


SUM5*SUM5+U1(N)*U2(N)/LAM(N) 

00004400 

ISN 

0045 


SUM6'SUM6+U2(N)**2/LAM(N) 

00004500 

ISN 

0046 


WRITE (6 , 2 1 )N . SUM 1 , SUM2 , SUM3 . SUM4 , SUM5 , SUM6 

00004600 

ISN 

0047 

21 

FORMAT! 1H0.5X, 12. 6(G13 5.2X)) 

00004700 

ISN 

0048 

22 

CONTINUE 

00004800 

ISN 

0049 

23 

FORMAT! 1H0.42X, 'PARTIAL SUMS' , //////) 

00004810 



C 


00004900 



C 

EXACT SUMS 

00005000 



C 


00005100 

ISN 

0050 


SUM 1 * 1 +MSTAR 

00005200 

ISN 

0051 


SUM2 * 1 /3 +MSTAR+USTAR+2 *MSTAR*CSTAR 

00005300 

ISN 

0052 


SUMS'. 5+MSTAR*! 1 +CSTAR ) 

00005400 

ISN 

0053 


SUM4= 1 

00005500 

ISN 

0054 


SUMS' 5+CSTAR 

00005600 

ISN 

0055 


SUMS' 1 /3.+CSTAR*(1.+CSTAR) 

00005700 

ISN 

0056 


WRI TE ( 6 . 30 ) SUM 1 . SUM2 . SUMS , SUM 4 , SUM5 . SUM6 

00005800 

ISN 

0057 

30 

FORMAT! 1H0, //////, 8X,6(G 13 5,2X)) 

00005900 

ISN 

0058 


STOP 

00006000 

ISN 

0059 


END 

00006100 


ISN 0002 


FUNCTION CHARDT! ALF . MSTAR , dSTAR , CSTAR ) 

00006200 


C 


00006210 


C 

THIS FUNCTION COMPUTES THE EQUIVALENT OF EQ 2-19 

00006220 


C 


00006230 

ISN 0003 


IMPLICIT REAL*16(A-H.0-Z) 

00006300 

ISN 0004 


REAL* 16 MSTAR, USTAR 

00006400 

ISN 0005 


S'OSIN(ALF) 

00006500 

ISN 0006 


C'OCOS(ALF) 

00006600 

ISN 0007 


SH'OSINH(ALF) 

00006700 

ISN 0008 


CH'QCOSH(ALF) 

00006800 

ISN 0009 


C 1 = MST AR * ALF 

00006900 

ISN 0010 


C2'C1*CSTAR*ALF 

00007000 

ISN 0011 


C3=USTAR*ALF * *3 

00007100 

ISN 0012 


A11=C2*(C-CH)-C3*!S+SH)+C+CH 

00007200 

ISN 0013 


A12=C2*(S-SH)+C3*!C-CH)+S+SH 

00007300 

ISN 0014 


A21=C1*(CH-C)+C2*(S+SH)-S*SH 

00007400 

ISN 0015 


A22=C1*(SH-S)+C2*(CH-C)+C+CH 

00007500 

ISN 0016 


CHARDT' A 1 1*A22-A12*A21 

00007600 

ISN 0017 


RETURN 

00007700 

ISN 0018 


END 

00007800 
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ISN 0002 


SUBROUTINE E IGEN(NF , MSTAR . JSTAR . CSTAR . LAM, U 1 . U2 . U3 . U4 ) 

0000790C 


C 


00007910 


C 

THIS SUBROUTINE COMPUTES "NF" EIGENVALUES (ROOTS OF EO 2-19) AND 

00007920 


C 

MODAL PARAMETERS (EO'S 3-25) OF A CLAMPED BEAM WITH TIP BODY 

00007930 


C 


00007950 

ISN 0003 


IMPLICIT REAL*16(A-H,0-Z) 

00008000 

ISN 0004 


REAL* 16 MSTAR, JSTAR, LAM(NF), LEFT 

00008100 

ISN 0005 


DIMENSION U1(NF),U2(NF),U3(NF),U4(NF) 

00008200 


C 


00008320 


C 

ESTIMATE ROOT OF CHARACTERISTIC EOUATION (EO 2-19) 

00008330 


C 


00008340 


C 

N IS THE NUMBER OF ROOTS FOUND 

00008400 


C 


00008500 

ISN 0006 


N=0 

00008600 

ISN 0007 


5TEP= 1 00-02 

00008700 

ISN 0003 


LEFT =0 000 

00008800 

ISN 0009 

1 

RIGHT =LEFT+STEP 

00008900 


C 


00009000 


C 

COMPARE THE SIGNS OF THE FUNCTION "CHARDT" AT THE TWO POINTS 

00009100 


C 


00009200 

ISN 0010 


IF(CHARDT( LEFT, MSTAR, JSTAR .CSTAR )*CHARDT (RIGHT .MSTAR, JSTAR, CSTAR) 

00009300 



1 LE 0 000) GO TO 2 

00009400 

ISN 0012 


LEFT = RIGHT 

00009500 

ISN 0013 


GO TO 1 

00009600 

ISN 0014 

2 

N=N+1 

00009700 

ISN 0015 


ALF=RIGHT 

00009800 


C 


00009900 


C 

IMPROVE ROOT ESTIMATE WITH NEWTON-RAPHSON ITERATION 

00010000 


C 


00010100 

ISN 0016 


NIT=0 

00010200 

ISN 0017 

3 

S=0SIN( ALF ) 

00010300 

ISN 0018 


C=OCOS( ALF ) 

00010400 

ISN 0019 


SH=0SINH(ALF) 

00010500 

ISN 0020 


CH=QCOSH( ALF ) 

00010600 

ISN 0021 


C 1 'MSTAR* ALF 

00010700 

ISN 0022 


C2=C1*CSTAR*ALF 

00010800 

ISN 0023 


C3=JSTAR*ALF**3 

00010900 

ISN 0024 


A11=C2*(C-CH)-C3*(S+SH)+C+CH 

00011000 

ISN 0025 


A 12=C2* ( S"SH)+C3* ( C-CH ) +S+SH 

00011100 

ISN 0026 


A2 1 =C1 * (CH-C )+C2*( S+SH) -S+SH 

00011200 

ISN 0027 


A22=C1*(SH-S)+C2*(CH-C)+C+CH 

0001 1300 

ISN 0028 


A11P=MSTAR*CSTAR*(2 000* ALF * (C-CH ) - ALF * *2 * ( S+SH) ) -C3* (C+CH ) 

00011400 



1 -3 000 *JS TAR* ALF* *2* (S+SH) -S+SH 

0001 1500 

ISN 0029 


A12P=MSTAR*CSTAR*(2 000*ALF* (S-SH)+ALF**2* (C-CH) )-C3*( S+SH) 

00011600 



1 +3 OQO*JSTAR*ALF* *2*(C-CH)+C + CH 

00011700 

ISN 0030 


A21P=C1*(S+SH)+C2*(C+CH)+2 OQO*MSTAR 

0001 1800 



1 *CSTAR * ALF* (S+SH) -C+CH+MSTAR* (CH-C) 

00011900 

ISN 0031 


A22P=MSTAR* ( SH- S )+C 1 * (CH-C )+C2* ( SH+S ) 

00012000 



1 +2 000*MSTAR ‘CSTAR* ALF* (CH-C) -S+SH 

00012100 

ISN 0032 


F=A11*A22-A12*A21 

00012200 

ISN 0033 


FP=A1 1P*A22+A1 1*A22P-A12P*A21-A12*A21P 

00012300 

ISN 0034 


DA-F/FP 

00012400 

ISN 0035 


ALF=ALF -DA 

00012500 
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ISN 

0036 


IF(QABS(DA/ALF) LE 1.00-09) GO TO 6 

00012600 

ISN 

0038 


NIT = NIT+ 1 

00012700 

ISN 

0039 


IF(NIT GT. 10) GO TO 4 

00012800 

ISN 

0041 


GO TO 3 

00012900 

ISN 

0042 

4 

WRITE(6.5) NIT.N.OA.ALF 

00013000 

ISN 

0043 

5 

FORMAT! 1H0.5X. 'NEWTON RAPHSON ITERATION FAILED TO CONVERGE 

IN' . IX, 00013100 




1 12,' ITERATIONS. ROOT # ' . I2.2X, 'DA«' , E 13 5,2X. ' ALF* ' , E 13 

5) 00013200 

ISN 

0044 

6 

S*QSIN( ALF ) 

00013300 

ISN 

0045 


C=QCOS( ALF ) 

00013400 

ISN 

0046 


SH=QSINH(ALF ) 

00013500 

ISN 

0047 


CH=QCOSH(ALF) 

00013600 

ISN 

0048 


C1=MSTAR*ALF 

00013700 

ISN 

0049 


C2=C1*CSTAR*ALF 

00013800 

ISN 

0050 


C3=JSTAR*ALF‘«3 

00013900 

ISN 

0051 


A11-C2*(C-CH)-C3*( S+SH )+C+CH 

00014000 

ISN 

0052 


A12-C2*(S-SH)+C3*(C-CH)+S+SH 

00014100 

ISN 

0053 


LAM(N)=ALF**4 

00014200 

ISN 

0054 


BETA = -A1 1/A 12 

00014300 

ISN 

0055 


T 1=ALF*(SH+S+BETA*(CH-C) ) 

00014400 

ISN 

0056 


T2=CH-C+BETA*(SH-S) 

00014500 

ISN 

0057 


T3 = ( SH-S + BETA * (CH+C-2 .000) J/ALF 

00014600 

ISN 

0058 


T4=(SH-S+BETA‘(CH+C))/ALF-(CH+C-2 000+BETA*(S+SH) )/ALF‘*2 

00014700 



C 


00014800 



C 

CALCULATE MODAL PARAMETERS WITH NON NORMALIZED EIGENFUNCTIONS 

00014900 



C 


00015000 

ISN 

0059 


U1(N)-T1 

00015100 

ISN 

0060 


U2(N) = T2+CSTAR*T 1 

00015200 

ISN 

0061 


U3(N)-T3+MSTAR*T2+MSTAR*CSTAR*T 1 

00015300 

ISN 

0062 


U4(N)=T4+MSTAR*(1 OOO+CSTAR )*T2+(MSTAR*CSTAR+USTAR ) *T 1 

00015400 



C 


00015500 



C 

CALCULATE INTEGRAL OF SQUARE OF EIGENFUNCTION 

00015600 



C 


00015700 

ISN 

0063 


C1 = 1 OOO+BETA* *2 

00015800 

ISN 

0064 


C2=BETA* *2- 1 000 

00015900 

ISN 

0065 


T5=(C1«CH+2 000*BETA*SH)*SH/(2 000*ALF) 

00016000 

ISN 

0066 


T6-C2«SH*C/ALF-C1*CH‘S/ALF-2.000*BETA*SH*S/ALF 

00016100 

ISN 

0067 


T7=S* (2 000*BETA*S-C2*C)/(2.000‘ALF)+1 OQO 

00016200 

ISN 

0068 


VALINT=T5+T6+T7 

00016300 



C 


00016400 



C 

CALCULATE NORM OF EIGENFUNCTION & MODAL PARAMETERS FOR NORMALIZED 00016500 



C 

EIGENFUNCTIONS 

00016600 



C 


00016610 

ISN 

0069 


P2=VALINT+MSTAR*T2**2*JSTAR*T1»*2+2 000*MSTAR*CSTAR*T 1 *T2 

00016700 

ISN 

0070 


P'0S0RT(P2) 

00016800 

ISN 

0071 


U1 (N) -U1 (N)/P 

00016900 

ISN 

0072 


U2(N)=U2(N)/P 

00017000 

ISN 

0073 


U3(N)*U3(N)/P 

00017100 

ISN 

0074 


U4(N)-U4(N)/P 

00017200 

ISN 

0075 


IF(N .EO NF) RETURN 

00017300 

ISN 

0077 


LEFT=RIGHT 

00017400 

ISN 

0078 


GO TO 1 

00017500 

ISN 

0079 


END 

00017600 
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OUTPUT FROM "EIGENVALUE" -MODAL PARAMETER PROGRAM 


DATA FROM NAMELIST INPUT 

MSTAR- 2 OOOO USTAR* 0 0280 CSTAR* 0 1000 
NF* 10 


■EIGENVALUES" & MODAL PARAMETERS 


N 

LAM 

U1 

U2 

U3 

U4 

1 

1 0310 

0 9087 

0 6760 

1 

56911 

1 6540 

2 

143 31 

-4 8354 

-0 1266 

0 

52240 

14854 

3 

1220 0 

5 0703 

-0 0027 

0 

29800 

505870-01 

4 

5231 5 

-4 9666 

0 0552 

0 

22042 

259090-01 

5 

16775 

3 5599 

-0 0608 

0 

17072 

150010-01 

S 

42936 

-2 6385 

0 0551 

0 

13693 

951230-02 

7 

93095 

2 0584 

-0 0485 

0 

11354 

650020-02 

8 

178940+06 

-1 6739 

0 0427 

0 

09673 

470290-02 

9 

314510+06 

1 4044 

-0 0380 

0 

08415 

355330-02 

10 

516170+06 

-1 2066 

0 0341 

0 

07442 

277650-02 



N 

N 

MODAL PARAMETER 
PARTIAL 

N 

IDENTITIES 

SUMS 

N „ 2 

N „ i i 

N ,, 2 


£ 4, 

k-1 

£ -L 

k-1 

/ u«, u.. 

V 1 Ik 

^ u 1k u 2k 

U 2k 

N 

3k 4k 
k-1 

■ * X 

k-1 *k 

h H 

1 

2 4621 

2.7356 

2 5953 

80094 

59581 

44322 

2 

2 7350 

2 7577 

2 6729 

96409 

60008 

44333 

3 

2 8238 

2 7602 

2 6879 

99429 

60007 

44333 

4 

2 8724 

2 7609 

2 6936 

99901 

60002 

44333 

5 

2 9015 

2 7611 

2 6962 

99976 

60001 

44333 

6 

2 9203 

2.7612 

2 6975 

99993 

60000 

44333 

7 

2 9332 

2 7613 

2 6983 

99997 

60000 

44333 

a 

2 9425 

2 7613 

2 6987 

99999 

60000 

44333 

9 

2 9496 

2.7613 

2 6990 

99999 

60000 

44333 

10 

2 9552 

2 7613 

2 6992 

1 0000 

60000 

44333 


1 +m* 

1/3 + m* 

+ J* + 2m*c* 

ASYMPTOTIC VALUES 
1/2 + m* +m*c* 1 0 

1/2 + c* 

1/3 + c* +c* : 


3 0000 

2 7613 

2 7000 

1 0000 

60000 

44333 


80 



APPENDIX C 


PLANAR DYNAMICS FORTRAN PROGRAM LISTING 


ISN 0002 
ISN 0003 
ISN 0004 
ISN 0005 


ISN 0006 
ISN 0007 
ISN 0008 
ISN 0009 
ISN 0010 


C********************************************************************** 

C* THIS PROGRAM SOLVES THE EQUATIONS OF MOTION FOR THE PLANAR DYNAMICS 
C* OF A FLEXIBLE BEAM WITH RIGID BODIES ATTACHED TO EACH END THE 
C* BEAM DEFORMATION IS EXPANDED IN TERMS OF THE MODES APPROPRIATE TO 
C* A CLAMPED BEAM WITH TIP BODY ONLY THE ATTITUDE AND BENDING COOR- 
C* DINATES ARE CALCULATED ( WRITTEN BY JOEL STORCH & STEPHEN GATES 
C* BASED UPON C S 0 L REPORT # R-1629 MAY 1983 ) 

C******* *************************************************************** 

c 

C NOTE ARRAYS DIMENSIONED FOR A MAXIMUM OF 20 CLAMPED-TIP BODY MODES 
C 

IMPLICIT REAL’8(A-H.0-Z) 

REAL* 16 LAM , U1 ,U2 , U3 , U4 , MSTAR , JSTAR , CST 
REAL *8 M,L,MO, IO.MT.IT, Ml .MU0.MU1, JO.NA 

DIMENSION LAM(20).U1(20).U2(20).U3(20).U4(20),A1V(20),AS(231 ), 

1 B5(231),Z(21,21),EV(21) , WK( 483),FV(21),CN5(20), CN6( 20) , CN7 ( 20) , 

2 CN8(20),F0(2),FA(21),NA(21),W1(21),W2(21).Y(42).YD0T(42) 

COMMON /DIM/ NF.NFP1.NSZ 

COMMON /EXTFC/ CN 1 , CN2 , CN3 , CN4 , CN5 , CN6 , CN7 
COMMON /NLKTC/ CN9.U3.CN8 
COMMON /STATE/ X ( 2 1 ) . XD0T(2 1 ) 

NAMELIST /INPUT/ MO, 10. A 1 , A2. L. PHO , C . MT , IT . NF , El . THETA , THETAD , 

1 DT.TSTOP.TPRT 






INPUT - OUTPUT FILES 
FILE #5 NAMELIST INPUT FILE 
FILE H 6 PRINTED OUTPUT FILE 

FILE #8 INPUT FILE OF EIGENVALUES & MODAL PARAMETERS FOR BEAM WITH 
TIP BODY 

****** ***+*******DESCRIPTION OF /INPUT/ LIST ITEMS***************** 


"MO” 

" 10 " 

“ A 1 " , 

"RHO" 

"El" 

»C" 

"MT" 

“IT" 

"NF " 


ORBITER MASS 

ORBITER MOMENT OF INERTIA ABOUT ITS MASS CENTER 
PERPENDICULAR TO PLANE OF MOTION 

A2 " VECTOR FROM ORBITER MASS CENTER TO BEAM ATTACHMENT POINT 
(EXPRESSED IN ORBITER FRAME) 

BEAM LENGTH 

MASS PER UNIT LENGTH OF BEAM 
BENOING STIFFNESS IN BEAM 
TIP BODY MASS CENTER OFFSET 
TIP BODY MASS 

TIP BODY MOMENT OF INERTIA ABOUT ITS MASS CENTER FOR AXIS 
PERPENDICULAR TO PLANE OF MOTION 

NUMBER OF CANTILEVERED-TIP BODY MODES TO BE RETAINED 
FOR EXPANSION OF BEAM ELASTIC DISPLACEMENT 
"THETA" ORBITER INITIAL ATTITUDE (DEG ) 

"THETAD" ORBITER INITIAL ATTITUDE RATE (DEG /SEC) 

"DT" INTEGRATION TIME STEP 
"TSTOP" FINAL TIME 

"TPRT " TIME INTERVAL BETWEEN CONSECUTIVE PRINTOUTS 
***************************************************************** 


00000100 

*00000200 

*00000300 

*00000400 

*00000500 

*00000600 

*00000700 

*00000800 

*00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

00002100 

00002200 

00002300 

00002400 

00002500 

*00002600 

00002700 

00002800 

00002900 

00003000 

00003100 

00003200 

00003300 

00003400 

00003500 

00003600 

00003700 

00003800 

00003900 

00004000 

00004100 

00004200 

00004300 

00004400 

00004500 

00004600 

00004700 

00004800 

00004900 

00005000 

00005100 

00005200 

00005300 


81 




C 




00005400 

ISN 0011 


READ(5, INPUT) 



00005500 

ISN 0012 


WRITE(6,7) 



00005700 

ISN 0013 


WRITE(6, 1) MO, 10 



00005800 

ISN 0014 

1 

FORMAT( 1H0, 'MO = ' . 1PE 17 .8, 5X. ' 10 = 

'. 1PE17.8) 


00005900 

ISN 0015 


WRITE (6 , 2 ) MT , IT 



00006000 

ISN 001S 

2 

F0RMAT(1H , 'MT = ' , 1PE17 8,5X,'IT * 

' , 1PE17.8) 


00006 100 

ISN 0017 


WRITE(6,3) A1.A2 



00006200 

ISN 0018 

3 

FORMAT ( 1H , ' A 1 = '.1PE17 8,5X,'A2 = 

' . 1PE17 8) 


00006300 

ISN 0019 


WRITE(6,4) L,C 



00006400 

ISN 0020 

4 

FORMAT ( 1H ,'L - ' , 1PE 17 8,5X,'C - 

', 1PE17 8) 


00006500 

ISN 0021 


WRITE(6,5) RHO.EI 



00006600 

ISN 0022 

5 

F0RMAT(1H ,'RHO ='.1PE17 8,5X.'EI = 

' . 1PE17 8) 


00006700 

ISN 0023 


WRITE{6,6) DT.TPRT 



00006710 

ISN 0024 

6 

FORMAT ( 1H ,'DT = '.1PE17 8,5X.'TPRT 

»' . 1PE16 8) 


00006720 

ISN 0025 


WRITE(6,8) THETA. THETAO 



00006800 

ISN 0026 

7 

F0RMAT(1H1,13X, , DATA FROM NAMELIST INPUT') 


00007000 

ISN 0027 

8 

FORMAT ( 1H , 'THETA °',1PE15 8 , 3X , ' THETADOT *'1PE15 8) 


00007010 


C 




00007100 


C 

TEST FOR DATA CONSISTENCY WITH VALUES 

ON DISC DATA SET 


00007200 


C 




00007300 

ISN 0028 


READ(8) NFT.MSTAR.USTAR.CST 



00007400 

ISN 0029 


D 1 =DABS(DBLEO(MSTAR ) -MT/(RHO*L ) ) 



00007500 

ISN 0030 


02>DABS(0BLEQ(JSTAR)-!IT+MT*C**2)/(RH0*L**3)) 


00007600 

ISN 0031 


D3=DABS(DBLE0(CST)-C/L) 



00007700 

ISN 0032 


D 1 = D 1/DBLE0( MSTAR ) 



00007800 

ISN 0033 


D2-02/DBLEQ(USTAR) 



00007900 

ISN 0034 


03*03/DBLEQ(CST) 



00008000 

ISN 0035 


IER-0 



00008100 

ISN 0036 


I F (NFT GE NF) GO TO 11 



00008200 

ISN 0038 


IER= 1 



00008300 

ISN 0039 


WRITE(G.IO) NF.NFT 



00008400 

ISN 0040 

10 

FORMAT! 1H0.5X.'* • * FATAL ERROR ♦ ♦ 

* ' , 2X . 13 . ' CANTILEVERED' 


00008500 



1 ' TIP-BODY MODES REQUESTED DATA ON 

DISC ONLY FOR', 13. 


00008600 



2 ' MODES') 



00008700 

ISN 0041 

11 

I F( D 1 LE 00100) GO TO 13 



00008800 

ISN 0043 


IER»1 



00008900 

ISN 0044 


WRITE (6 . 12) MSTAR 



00009000 

ISN 0045 

12 

FORMAT! 1H0.5X, '* ♦ ♦ FATAL ERROR * * 

• ' ,2X, 'MSTAR(DISC)*' .G13 

5) 

00009100 

ISN 0046 

13 

IF!02 LE 00100) GO TO 15 



00009200 

ISN 0048 


IER* 1 



00009300 

ISN 0049 


WRITE ( 6 , 14) JSTAR 



00009400 

ISN 0050 

14 

FORMAT! 1H0.5X, '* * * FATAL ERROR* * 

* ' ,2X,'dSTAR!DISC)='.G13 

5) 

00009500 

ISN 0051 

15 

IF(D3 .LE 00100) GO TO 17 



00009600 

ISN 0053 


IER* 1 



00009700 

ISN 0054 


WRITE(6, 16) CST 



00009800 

ISN 0055 

16 

FORMAT! 1H0.5X, '* * ♦ FATAL ERROR* * 

* '.2X,'CSTAR(DISC)*',G13 

5) 

00009900 

ISN 0056 

17 

IF( IER NE 0) STOP 



00010000 


C 




00010100 


C 

COMPUTE CONSTANT QUANTITIES 



00010200 


C 




00010300 

ISN 0058 


R5=EI/(RH0*L**4) 



00010400 

ISN 0059 


FC« 1591549*DS0RT(R5) 



00010500 

ISN 0060 


M1=RH0*L+MT 



00010600 

ISN 0061 


M*M0+M1 



00010700 

ISN 0062 


MU0=M0/M 



00010800 

ISN 0063 


MU1 =M1/M 



00010900 

ISN 0064 


B1=! 5*RH0*L**2+MT*(L+C))/M1 



00011000 


82 



ISN 

0065 


d0=RH0*L**3/3 OOO^IT+MT* (C+L ) * *2 

00011100 

ISM 

0066 


R 1=RH0*L/M 

0001 1200 

ISN 

0067 


R2=A 1/L 

00011300 

ISN 

0068 


R3=A2/L 

0001 1400 

ISN 

0069 


R4=B1/L 

00011500 

ISN 

0070 


CSTAR=C/L 

00011600 

ISN 

0071 


AOO=(IO+dO)/(RHO*L**3)+MUO*MU1*(R2**2+R3**2+2 0D0*R2*R4 )/R 1 

00011700 




1 -MU1 **2*R4**2/R1 

00011800 



C 


00011900 



C 

READ & PRINT EIGENVALUES ANO MODAL PARAMETERS FOR BEAM+TIP BODY 

00012000 



C 


00012100 

ISN 

0072 


WRITE (6 , 23 ) 

00012110 

ISN 

0073 


WRITE(6,21) 

00012200 

ISN 

0074 


DO 20 N=1,NF 

00012300 

ISN 

0075 


READ(8 ) LAM(N),U1(N).U2(N),U3(N),U4(N) 

00012400 

ISN 

0076 

20 

WRITE (6. 22) N,LAM(N),U1(N),U2(N),U3(N),U4(N) 

00012500 

ISN 

0077 

21 

FORMAT ( 1H0.T4, 'N' ,T 12, 'LAM' . T27 , ' U 1' , T40, ' U2 ' , T49 , ' U3 ' , T66 , 'U4 ' ) 

00012600 

ISN 

0078 

22 

FORMAT ( 1H ,T3,I2,T8,G12 5.T24.F9 4.T36.F9 4.T47.F9 5.T60.G13 5) 

00012800 

ISN 

0079 

23 

FORMAT( 1H0,/, 11X, 'EIGENVALUES & MODAL PARAMETERS FROM DISC FILE') 

00012810 

ISN 

0080 


WRITE { 6 , 24 ) 

00012820 

ISN 

0081 

24 

FORMAT( 1H0./.7X, 'SYSTEM NATURAL FREQUENCIES') 

00012830 



C 


00012900 



C 

COMPUTE TERMS All 

00013000 



C 


00013100 

ISN 

0082 


DO 30 K= 1 , NF 

00013200 

ISN 

0083 

30 

A 1 V(K ) =MU0*R2*U3(K )+U4 (K ) -MU1 *R4*U3(K ) 

00013300 



C 


00013400 



C 

STORE "A” IN SYMMETRIC STORAGE MODE - "AS" 

00013500 



C 

STORF "K" IN SYMMETRIC STORAGE MODE - "BS" 

00013600 



C 


00013700 

ISN 

0084 


AS( 1 )-A00 

00013800 

ISN 

0085 


N= 1 

00013900 

ISN 

0086 


DO 41 1=1, NF 

00014000 

ISN 

0087 


N=N-M 

00014100 

ISN 

0088 


AS(N)=A1V(I) 

00014200 

ISN 

0089 


DO 40 J»1,I 

00014300 

ISN 

0090 


N=N+1 

00014400 

ISN 

0091 


AS(N)=-R1*U3(d)*U3(I) 

00014500 

ISN 

0092 


IF( I EO d) AS(N)=AS(N)+1 ODO 

00014600 

ISN 

0094 

40 

CONTINUE 

00014700 

ISN 

0095 

41 

CONTINUE 

00014800 

ISN 

0096 


BS( 1 )=0 ODO 

00014900 

ISN 

0097 


N= 1 

00015000 

ISN 

0098 


DO 51 1= 1 ,NF 

00015100 

ISN 

0099 


N=N+ 1 

00015200 

ISN 

0100 


BS(N)=0 000 

00015300 

ISN 

0101 


DO 50 d=1,I 

00015400 

ISN 

0102 


N=N+1 

00015500 

ISN 

0103 


BS(N) =0 ODO 

00015600 

ISN 

0104 


I F ( I EO d) BS(N) =LAM( I ) 

00015700 

ISN 

0106 

50 

CONTINUE 

00015800 

ISN 

0107 

51 

CONTINUE 

00015900 



C 


00016000 



C 

GET EIGENVALUES ANO EIGENVECTORS FOR 2,3, .NFP1 DEGREES OF FREEDOMOOO 16 100 



C 


00016200 

ISN 

0108 


NFP 1 =NF+ 1 

00016300 

ISN 

0109 


NSZ=2»NFP1 

00016400 

ISN 

0110 


DO 66 N=2 , NFP 1 

00016500 


83 



ISN 

0111 


CALL EIGZS(BS.AS.N, 1.EV.Z.21.WK.IER) 

00016600 

ISN 

0112 


IF( IER EQ. 0) GO TO 62 



00016700 

ISN 

0114 


WRITE(6.61) N.IER 



00016800 

ISN 

01 15 

61 

FORMAT* 1H0.5X, 'ERROR IN EIGENVALUE 

EXTRACTION MATRIX ORDER"' 

.12. 00016900 




1 3X. 'IER=' .13) 



00017000 

ISN 

0116 


STOP 



00017100 

ISN 

0117 

62 

WRITE ( 6 , 63 ) N 



00017200 

ISN 

0118 

63 

FORMAT ( 1H0, 12 , ' OEGREES OF 

FREEDOM 

IN EIGENVALUE PROBLEM') 

00017300 

ISN 

0119 


DO 64 1=1, N 



00017500 

ISN 

0120 


WS=FC«OSORT(EV(I)) 



0001 7600 

ISN 

0121 

64 

WRITE ( 6 . 65 ) I.WS 



00017700 

ISN 

0122 

65 

F0RMAT(1H , 2X . ' MOOE ' , 12 , 2X , 

'SYSTEM 

FREQUENCY(HZ )='.G13 5) 

00017800 

ISN 

0123 

66 

CONTINUE 



00017900 



C 




00018000 



C 

NORMALIZE EIGENVECTORS 



00018100 



C 




00018200 

ISN 

0124 


DO 70 N= 1 ,NFP1 



00018300 

ISN 

0125 

70 

CALL N0RM(AS.NFP1,Z(1,N)) 



00018400 



C 




00018500 



C 

CALCULATE CONSTANTS IN MOTION 

EQUATIONS 

00018600 



C 




00018700 

ISN 

0126 


DO 80 1=1, NFP 1 



00018800 

ISN 

0127 

80 

FV(I)=R5 , EV(I) 



00018900 

ISN 

0128 


CN 1 = RHO* L* *3 



00019000 

ISN 

0129 


CN2=MU 1 * ( A 1+B 1 )/CN1 



00019100 

ISN 

0130 


CN3=MU 1 * A2/CN 1 



00019200 

ISN 

0131 


CN4=( 1 0D0+CSTAR+MU0‘R2-MU1 

♦R4)/(RH0*L‘*2) 

00019300 

ISN 

0132 


CN9= 1 ODO/(M‘L) 



00019400 

ISN 

0133 


DO 90 I = 1 , NF 



00019500 

ISN 

0134 


CN5( I )=U1(I)/CN1 



00019600 

ISN 

0135 


CN6( I )=U2( I )/(RH0*L**2 ) 



00019700 

ISN 

0136 


CN7 ( I ) =U3( I ) *CN9 



00019800 

ISN 

0137 

90 

CN8( I ) =MU0*R3*U3( I ) 



00019900 



C 




00020000 



c 

SET INITIAL CONDITIONS 



00020100 



c 




00020200 

ISN 

0138 


T =0 ODO 



00020300 

ISN 

0139 


CPRT =0 ODO 



00020400 

ISN 

0140 


X( 1 ) -THETA * 0174532 



00020500 

ISN 

0141 


XDOT( 1 )=THETAD* 0174532 



00020600 



c 




00020700 



c 

INITIAL DEFORMATION AND RATE 

ARE SET 

TO ZERO 

00020800 



c 




00020900 

ISN 

0142 


DO 100 1=1. NF 



0002 1000 

ISN 

0143 


X ( I 1 )=0.0D0 



00021100 

ISN 

0144 

100 

XDOT( 1+ 1 ) =0.000 



00021200 



c 




00021300 



c 

CALCULATE Y AT T=0 



00021400 



c 




00021500 

ISN 

0145 


DO 112 1=1, NFP 1 



00021600 

ISN 

0146 


II = I 



00021700 

ISN 

0147 


W1(I)«0 ODO 



00021800 

ISN 

0148 


W2( I )=0 000 



00021900 

ISN 

0149 


DO 111 J= 1 , NFP 1 



00022000 

ISN 

0150 


JvJ=vJ 



00022100 

ISN 

0151 


I F ( I GE. 0) GO TO 110 



00022200 

ISN 

0153 


II=U 



00022300 

ISN 

0154 


JJ=I 



00022400 


84 



ISN 0155 

1 10 

LL = II + ( II - 1 )/2+UJ 

00022500 

ISN 0156 


W1(I)=W1(I)+AS(LL)’X(U) 

00022600 

ISN 0157 

111 

W2( I ) =W2( I)+AS(LL ) *XDOT( J ) 

00022700 

ISN 0158 

1 12 

CONTINUE 

00022800 

ISN 0159 


DO 114 1=1, NFP 1 

00022900 

ISN 0160 


Y ( I ) =0 ODO 

00023000 

ISN 0161 


Y ( NFP 1+1 ) =0 ODO 

00023100 

ISN 0162 


DO 113 J= 1 , NFP 1 

00023200 

ISN 0163 


Y(I)=Y(I)+Z(J,I)*W1(J) 

00023300 

ISN 0164 

113 

Y(NFP1+I)=Y{NFP1+I)+Z(J,I)*W2(U) 

00023400 

ISN 0165 

114 

CONTINUE 

00023500 

ISN 0166 


WRITE(6, 115) 

00023510 

ISN 0167 

115 

FORMAT ( 1H0,/,20X, 'TIME RESPONSE' ) 

00023520 


C 


00023600 


C 

CALCULATE EXTERNAL FORCES ON ORBITER AND TIP BODY 

00023700 


C 


00023800 

ISN 0168 

120 

CALL ORBFOR(T, FO.GO) 

00023900 

ISN 0169 


CALL TPBFOR(T,FP,GP) 

00024000 

ISN 0170 


CALL EXTF(FO,GO,FP,GP,FA) 

00024100 


C 


00024200 


C 

CALCULATE NON-LINEAR TERMS 

00024300 


C 


00024400 

ISN 0171 


CALL NLKT(FO.NA) 

00024500 

ISN 0172 


00 131 1=1, NFP 1 

00024600 

ISN 0173 


W1 ( I ) =0 ODO 

00024700 

ISN 0174 


DO 130 J= 1 ,NFP 1 

00024800 

ISN 0175 

130 

W1( I )=W1(I ) + Z(U, I )*(FA(J)+NA(U) ) 

00024900 

ISN 0176 

131 

CONTINUE 

00025000 


C 


00025100 


C 

CALCULATE “YDOT" 

00025200 


C 


00025300 

ISN 0177 


DO 140 1=1, NFP 1 

00025400 

ISN 0178 


YDOT(I )=Y(NFP1+I ) 

00025500 

ISN 0179 


YD0T(NFP1+I )=W1(I)-FV(I)*Y(I) 

00025600 

ISN 0180 

140 

CONTINUE 

00025700 


C 


00025800 


C 

INTEGRATE DIFFERENTIAL EQUATIONS IN FIRST ORDER FORM 

00025900 


C 

- 

00026000 

ISN 0181 


CALL ODESLV(NSZ,Y, YOOT.DT ) 

00026100 

ISN 0182 


T=T+DT 

00026200 

ISN 0183 


CPRT =CPRT +DT 

00026300 

ISN 0184 


DO 151 1=1, NFP 1 

00026400 

ISN 0185 


X(I)=0 ODO 

00026500 

ISN 0186 


XDOT ( I ) =0 ODO 

00026600 

ISN 0187 


DO 150 U= 1 . NFP 1 

00026700 

ISN 0188 


X(I)=X(I)+Z(I,J)*Y(J) 

00026800 

ISM 0189 

150 

XDOT( I )=XDOT (I)+Z(I,U)*Y(NFP1+<J) 

00026900 

ISN 0190 

151 

CONTINUE 

00027000 

ISN 0191 


THETA=X( 1 )*57 29578 

00027100 

ISN 0192 


THETAD = XDOT( 1 )*57 29578 

00027200 


C 


00027300 


C PRINT OUTPUT 

00027400 


C 


00027500 

ISN 0193 


IF(CPRT LT TPRT ) GO TO 170 

00027600 

ISN 0195 


CPRT=0 ODO 

00027700 

ISN 0196 


WRITE(6, 160) T 

00027800 

ISN 0197 

160 

FORMAT(1HO, 'TIME = ' , F7 3) 

00027900 

ISN 0198 


WRITE(6, 161) THETA 

00028000 


85 



ISN 0199 


WRITE(6, 165) THETAD 

00028010 

ISN 0200 

161 

FORMAT ( 1H , 'THETA =',1PE17 8,' DEG') 

00028100 

ISN 0201 


WRITE (6 , 162) 

00028200 

ISN 0202 

162 

F0RMAT(1H0,' MODAL COORDINAT ES ' . 10X , ' DERI V MODAL COORDINATES') 

00028300 

ISN 0203 


DO 163 I=2,NFP1 

00028400 

ISN 0204 


J=I- 1 

00028500 

ISN 0205 

163 

WRITE(6, 164) J, X( I ) , XDOT (I ) 

00028600 

ISN 0206 

164 

FORMAT ( 1H .I2.2X.1PE17 8.10X.1PE17 8) 

00028700 

ISN 0207 

165 

FORMAT ( 1H ,'THETADOT «',1PE17 8.' DEG/SEC') 

00028710 


C 


00028800 

ISN 0208 

170 

I F ( T GE TSTOP ) STOP 

00028900 

ISN 0210 


GO TO 120 

00029000 

ISN 0211 


END 

00029100 


XSN 0002 


ISN 0003 
ISN 0004 
ISN 0005 
ISN 0006 
ISN 0007 
ISN 0008 
ISN 0009 
ISN 0010 


SUBROUTINE ORBFOR(T , FO.GO) 

C 

C THIS SUBROUTINE CALCULATES THE EXTERNAL FORCE "FO" ON THE ORBITER 

C (IN ORBITER FRAME) AMD NET MOMENT “GO" PERPENDICULAR TO PLANE OF 

C MOTION AT TIME “T“ 

C 

IMPLICIT REAL*8(A-H.0-Z) 

DIMENSION F0( 2 ) 

F0( 1 ) -0 ODO 
F0(2 )=0 ODO 
G0=0 000 

10 FORMAT ( 1H . 3X . ' F0= ' , 2E 13 5, 3X , ' G0= ' , E 13 5) 

RETURN 

END 


00029200 

00029300 

00029400 

00029500 

00029600 

00029700 

00029800 

00029900 

00030000 

00030100 

00030200 

00030300 

00030400 

00030500 


ISN 0002 SUBROUTINE TPBFOR( T , F , G ) 00030600 

C 00030700 

C THIS SUBROUTINE CALCULATES THE EXTERNAL FORCE ”F“ ON THE TIP BODY 00030800 
C (ACTING TRANSVERSE TO BEAM NEUTRAL AXIS) AND MOMENT "G” 00030900 

C PERPENDICULAR TO PLANE OF MOTION AT TIME “T" 00031000 

C 00031100 

ISN 0003 F=0 ODO 00031200 

ISN 0004 G=0 ODO 00031300 

ISN 0005 RETURN 00031400 

ISN 0006 END 00031500 


ISN 0002 


ISN 0003 
ISN 0004 
ISN 0005 
ISN 0006 
ISN 0007 
ISN 0008 
ISN 0009 
ISN 0010 
ISN 0011 


C 

c 

c 

c 

c 


10 


SUBROUTINE EXTF( FO, GO, FP . GP , F ) 

THIS SUBROUTINE ASSEMBLES THE VECTOR OF EXTERNAL FORCES “F" GIVEN 
BY EO (4-31) THE ORBITER FORCE AND MOMENT "FO’V’GO” AS WELL AS THE 
FORCE AND MOMENT ON THE TIP BODY "FP'V'GP" ARE INPUT 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION F0( 2 ) , F ( NFP 1 ) 

COMMON /DIM/ NF.NFP1.NSZ 

COMMON /EXTFC/ CN1 , CN2 , CN3 , CN4 , CN5(20 ) , CN6( 20) , CN7 ( 20 ) 

F( 1 )=G0/CN1-CN2»F0(2)+CN3*F0( 1 )+GP/CN1+CN4*FP 
DO 10 1=1, NF 

F( 1+ 1 )=CN5( I ) *GP+CN6( I ) *FP -CN7 (I)*(F0(2)+FP) 

RETURN 

END 


00031600 
0003 1 700 
00031800 
00031900 
00032000 
00032100 
00032200 
00032300 
00032400 
00032500 
00032600 
00032700 
00032800 
00032900 
00033000 
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ISN 0002 


SUBROUTINE NLKT(FO.N) 

00033100 


C 


00033200 


C 

THIS SUBROUTINE ASSEMBLES THE VECTOR OF NON-LINEAR TERMS "N" 

00033300 


C 

GIVEN BY EQ (4-32) THE EXTERNAL FORCE ON THE ORBITER "FO" IS 

INPUT 00033400 


C 


00033500 

ISN 0003 


IMPLICIT REAL*8(A-H,0-Z) 

00033600 

ISN 0004 


REAL*8 N 

00033700 

ISN 0005 


REAL* 16 U3 

00033800 

ISN 000G 


COMMON /NLKTC/ CN9 . U3(20) , CN8(20) 

00033900 

ISN 0007 


COMMON /STATE/ X( 2 1 ) , XDOT( 2 1 ) 

00034000 

ISN 0008 


COMMON /DIM/ NF.NFP1.NSZ 

00034100 

ISN 0009 


DIMENSION F0( 2 ) , N(NFP 1 ) 

00034200 

ISN 0010 


N( 1 ) =0 ODO 

00034300 

ISN 0011 


DO 10 K=1,NF 

00034400 

ISN 0012 

10 

N( 1 ) =N( 1)+0BLEQ(U3(K))*X(K+1) 

00034500 

ISN 0013 


N(l)=N(1 )*CN9*F0( 1 ) 

00034600 

ISN 0014 


THD2=XD0T( 1 )**2 

00034700 

ISN 0015 


DO 20 K=1,NF 

00034800 

ISN 0016 

20 

N(K+1)=CN8(K)*THD2 

00034900 

ISN 0017 


RETURN 

00035000 

ISN 0018 


END 

00035100 


ISN 

0002 


SUBROUTINE NORM(AS,N,X) 


00035200 



C 



00035300 



C 

THIS SUBROUTINE NORMALIZES THE EIGENVECTOR “X" WITH 

RESPECT TO THE 

00035400 



C 

SYMMETRIC POSITIVE DEFINITE MATRIX 'A' OF ORDER "N" 

STORED IN 

00035500 



C 

SYMMETRIC STORAGE MODE AS THE VECTOR "AS" 


00035600 



C 



00035700 

ISN 

0003 


IMPLICIT REAL*8(A-H,0-Z) 


00035800 

ISN 

0004 


DIMENSION A S ( 1 ) ,X( 1) 


00035900 

ISN 

0005 


S=AS(1)*X(1)**2 


00036000 

ISN 

0006 


NC= 1 


00036100 

ISN 

0007 


DO 10 1=2, N 


00036200 

ISN 

0008 


I M 1 = I - 1 


00036300 

ISN 

0009 


DO 20 J=1,IM1 


00036400 

ISN 

0010 


NC=NC+1 


00036500 

ISN 

0011 

20 

S = S*2 ODO* AS(NC )*X(I)*X(J) 


00036600 

ISN 

0012 


NC=NC+ 1 


00036700 

ISN 

0013 


S=S+AS(NC)*X(I)**2 


00036800 

ISN 

0014 

10 

CONTINUE 


00036900 

ISN 

0015 


S=DSQRT( S ) 


00037000 

ISN 

0016 


DO 30 I = 1 , N 


00037100 

ISN 

0017 

30 

X(I)=X(I)/S 


00037200 

ISN 

0018 


RETURN 


00037300 

ISN 

0019 


END 


00037400 


ISN 0002 


SUBROUTINE ODESLV(N, Y . DERI V.H) 

00037500 

ISN 0003 


IMPLICIT REAL*8(A-H.0-Z) 

00037600 

ISN 0004 


DIMENSION DERIV(N) , Y(N) . DERIV0(42 ) , BD 1 ( 42 . 2 ) , BD2 ( 42 . 2 ) ,BD3(42) 

00037700 

ISN 0005 

C 

OATA INTF/1/.C1/0 0/.C2/0 0/.C3/0 / 

00037800 

00037900 


C 

THIS SUBROUTINE INTEGRATES THE FIRST ORDER SYSTEM OF ORDINARY 

00038000 


C 

DIFFERENTIAL EQUATIONS "DY/DT=DERIV" BY THE ADAMS METHOD 

00038100 


C 

USING THIRD ORDER DIFFERENCES 

00038200 


C 

N- SIZE OF SYSTEM 

00038300 


C 

Y- VECTOR OF INITIAL VALUES ON INPUT "Y" IS OVERWRITTEN 

00038400 


C 

WITH THE NEW SOLUTION 

00038500 


C 

H- STEP SIZE 

00038600 
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c 00038700 


ISN 

0006 


I F ( N LE 42) GO TO 10 



00038800 

ISN 

0008 


WRITE (6 , 12) N 



00038900 

ISN 

0009 


STOP 



00039000 

ISN 

0010 

12 

FORMAt( 1H0.5X, 'ERROR IN SUBROUTINE 

**ODESLV»* CALLED WITH STATE 

00039100 




1SIZE = ' , 13 , ' EXCEEDS DIMENSION SIZE OF ARRAYS') 

00039200 

ISN 

0011 

10 

GO T0( 1000. 2000. 3000, 4000), INTF 


00039300 



C 




00039400 



C 

FIRST CALL TO ROUTINE - EULER 

INTEGRATION 

00039500 



C 




00039600 

ISN 

0012 

1000 

DO 20 1=1. N 



00039700 

ISN 

0013 

20 

DERIVO(I)=DERIV(I) 



00039800 

ISN 

0014 


INTF=2 



00039900 

ISN 

0015 


GO TO 5000 



00040000 



C 




00040100 



C 

SECOND CALL TO ROUTINE - FIRST 

ORDER 

DIFFERENCES 

00040200 



C 




00040300 

ISN 

0016 

2000 

DO 30 I =1 ,N 



00040400 

ISN 

0017 


BD 1 ( I , 1)=DERIV(I)-DERIV0(I) 



00040500 

ISN 

0018 


BD 1 ( I , 2 ) =BD 1(1,1) 



00040600 

ISN 

0019 

30 

DERIVO(I)=DERIV(I) 



00040700 

ISN 

0020 


C1= 5 



00040800 

ISN 

0021 


INTF =3 



00040900 

ISN 

0022 


GO TO 5000 



00041000 



c 




00041100 



c 

THIRD CALL TO ROUTINE - SECOND 

ORDER 

FIFFERENCES 

00041200 



c 




00041300 

ISN 

0023 

3000 

DO 40 I = 1 , N 



00041400 

ISN 

0024 


BD1(I.2)=DERIV(I )-DERIVO(I) 



00041500 

ISN 

0025 


BD2(I,1)=BD1(I.2)-BD1(I.1) 



00041600 

ISN 

0026 


B02( 1,2) =BD2( I . 1 ) 



00041700 

ISN 

0027 


DERIVO( I )=DERIV( I ) 



00041800 

ISN 

0028 

40 

BD1(I,1)=B01(I.2) 



00041900 

ISN 

0029 


INTF=4 



00042000 

ISN 

0030 


C2=5 0/12 0 



00042100 

ISN 

0031 


GO TO 5000 



00042200 



C 




00042300 



C 

ADAMS METHOD WITH 3RD ORDER DIFFERENCES 

00042400 



c 




00042500 

ISN 

0032 

4000 

DO 50 I = 1 , N 



00042600 

ISN 

0033 


BD1(I,2)=DERIV(I)-DERIV0(I) 



00042700 

ISN 

0034 


BD2(I,2)=BD1(I,2)‘BD1(I, 1) 



00042800 

ISN 

0035 


BD3 ( I ) = BD2( 1 , 2 ) -BD2 ( 1 , 1) 



00042900 

ISN 

0036 


DERIVO(I )=DERIV(I ) 



00043000 

ISN 

0037 


BD1 ( I , 1 )=BD 1(1,2) 



00043100 

ISN 

0038 

50 

BD2 ( I , 1 )=BD2( 1,2) 



00043200 

ISN 

0039 


C3=3 0/8 0 



00043300 

ISN 

0040 


GO TO 5000 



00043400 



C 




00043500 



C 

UPDATE VECTOR 'Y' 



00043600 



c 




00043700 

ISN 

0041 

5000 

DO 60 1=1, N 



00043800 

ISN 

0042 

60 

Y(I)=Y(I)+H*(DERIV(I)+C1*BD1(I,2 )+C2*BD2( 1 , 2 )+C3*B03( I ) ) 

00043900 

ISN 

0043 


RETURN 



00044000 

ISN 

0044 


END 



00044100 
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EXAMPLE PROBLEM PARAMETERS 


m o 

— 

98739.5 kg 

I o 

= 

9769869.5 kg-m' 

m t 

= 

875.32 kg 

x t 

= 

1400.512 kg-m 2 

a l 

= 

2.0 m 

a 2 

= 

0.0 m 

l 

= 

20.0 m 

c 

= 

2.0 m 

P 

= 

21.883 kg/m 

El 

= 

353520.0 N-m 2 

m* 

= 

2.0 

J* 

= 

0.028 

c* 

= 

0.1 


NAMELIST INPUT DATA 


&INPUT M0=98739 5,10=9769859 5,A1=2 ,A2=0 ,L=20 ,RH0=21 883, EI=353520 ,C=2 . 
MT=875 32, IT= 1400 5 12 , NF=3 . THETA=0 ,THETAD=0 ,DT = 01,TSTOP=1 ,TPRT= 02, SEND 
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OUTPUT DATA FROM PLANAR DYNAMICS PROGRAM 


4 

Excitation: G q = 4 x 10 N for all t >_ 0 

All units are metric (MKS) 


DATA FROM NAMELIST INPUT 


MO = 

9 87395000D+04 

10 = 

9 76986950D+06 

MT = 

8 75320000D+02 

IT = 

1 400512000+03 

A 1 = 

2 OOOOOOOOD+OO 

A2 = 

0 0 

L = 

2 00000000D+01 

C = 

2 OOOOOOOOD+OO 

RHO = 

2 18830000D+01 

El * 

3 53520000D+05 

DT = 

1 00000000D-02 

TPRT = 

2 000000000-02 

THETA 

*00 

THETADOT 

= 00 


EIGENVALUES & MODAL PARAMETERS FROM DISC FILE 

N LAM U1 U2 U3 U4 

1 1 0310 0 9087 0 S760 1 56911 1 6540 

2 143 31 -4 8354 -0 1266 0 52240 14854 

3 1220 0 6 0703 -0 0027 0 29800 505870-01 


SYSTEM NATURAL FREQUENCIES 

2 DEGREES OF FREEDOM IN EIGENVALUE PROBLEM 

MODE 1 SYSTEM FREQUENCY(HZ )* 0 

MODE 2 SYSTEM FREQUENCY(HZ )= 53106D-01 

3 DEGREES OF FREEDOM IN EIGENVALUE PROBLEM 

MODE 1 SYSTEM FREQUENCY(HZ )= 0 

MODE 2 SYSTEM FREQUENCY(HZ )= 53106D-01 

MODE 3 SYSTEM FREQUENCY ( HZ )= 60600 

4 DEGREES OF FREEDOM IN EIGENVALUE PROBLEM 

MODE 1 SYSTEM FREOUENCY(HZ ) = 0 

MODE 2 SYSTEM FREQUENCY(HZ ) = 53106D-01 

MODE 3 SYSTEM FREQUENCY(HZ )= 60600 

MODE 4 SYSTEM FREQUENCY(HZ ) = 1 7669 


TIME RESPONSE 

TIME = 0 020 

THETA = 3 5 1839 150D-05 DEG 

THETADOT = 4 69118867D-03 DEG/SEC 

MODAL COORDINATES DERIV MODAL COORDINATES 

1 -1 111891280-06 -1 48252 1 7 ID-04 

2 -1 23251214D-07 -1 64334952D-05 

3 -4 93395608D-08 -6 57860810D-06 

TIME = 0 040 

THETA = 1 759 19099D-04 DEG 

THETADOT = 9 38231697D-03 DEG/SEC 

MODAL COORDINATES DERIV MODAL COORDINATES 

1 -5 55938 160D -06 -2 964948560-04 

2 -6 15468068D-07 -3 276707300-05 

3 -2 44025899D-07 -1 281842280-05 
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End of Document 



